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INTRODUCTION 


The  crash  victim  simulator  (CVS)  computer  program  developed  by  Calspan 
is  a general  three-dimensional  simulator  (Ref.  1).  The  user  must  decide  the 
level  of  detail  (i.e.,  the  level  of  sophistication)  that  is  required  to  adequately 
account  for  all  aspects  of  a particular  problem  that  significantly  affect  the 
responses  of  interest.  In  the  crash  victim  problem  the  principal  responses  of 
interest  are  usually  those  associated  with  the  production  of  injuries.  The 
presently  reported  study  does  not  address  the  problem  of  identifying  injury 
criteria.  Rather,  the  general  objective  of  the  research  was  to  develop,  evaluate 
and  recommend  means  and  procedures  for  evaluating  accuracies  and  sensitivities 
of  various  measures  of  crash  victim  response  to  the  choice  of  and  to  the  numerical 
values  of  various  parameters  that  are  used  in  crash  victim  simulations  with  the 
CVS  computer  program. 

In  particular,  many  of  the  studies  were  confined  to  the  modeling 
aspects  related  to  a Part  572  dummy.  Extensive  work  to  define  the  parameters 
of  a 15  segment  Part  572  dummy  model  was  accomplished  in  two  closely  related 
research  programs  reported  in  References  1 and  2.  The  baseline  models  used  in 
this  study  stem  from  the  data  obtained  in  those  programs.  To  facilitate 
understanding  of  the  analyses  and  information  presented  in  this  report,  the 
reader  should  be  familiar  with  the  CVS  program  and  the  results  of  the  comparison 
studies  described  in  Reference  1. 

This  research  work  vras  divided  into  eight  (8)  principal  tasks  as 

follows : 

1.  Plan  of  Work  and  Methodology 

2.  Check  the  CVS  and  Input 

3.  Comparative  Study  of  Levels  of  Sophistication 

4.  Response  Measurq  Approximating  Function  Generator  (RMAFG) 

5.  Study  of  Responsfe  Sensitivity  to  Model  Parameter  Variations 

6.  Reduction  of  th^  Number  of  Inertial  Related  Parmaters 
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7.  Numerical  Verification  of  the  RMAFG 

8.  The  Main  Questions. 

Comments  on  the  Plan  of  Work  and  Methodology,  Task  1,  are  included 
in  the  discussions  of  the  individual  tasks.  The  original  plan  proved  overly 
ambitious.  Many  of  the  tasks  were  much  more  difficult  and  required  much  more 
time  than  had  originally  been  proposed.  As  a result  the  work  was  only 
partially  completed. 

Task  2,  Check  the  CVS  and  Input,  is  discussed  in  Section  2.  This  was 
principally  concerned  with  the  joint  algorithms  in  the  program  since  noticeable 
"drift"  of  the  joint  constraints  was  reported  by  some  users. 

Task  3,  Comparative  Studies  of  Levels  of  Sophistication,  addresses 
some  of  the  modeling  capabilities  of  the  CVS  program.  This  task  is  discussed 
in  Sections  3 and  6.  Section  6 also  addresses  Task  8,  the  Main  Questions. 

These  questions  are  given  in  Table  1-1  and  are  intimately  related  to  the  levels 
of  sophistication  that  were  to  be  considered  which  are  given  in  Table  1-2. 

Results  of  the  Study  of  Response  Sensitivity  to  Model  Parameter 
Variations  (Task  5)  are  presented  in  Section  3 where  a detailed  model  of  a 
dummy  shoulder  is  studied  and  in  Section  4 where  the  results  of  neck  and 
shoulder  models  are  used  to  demonstrate  the  Response  Measure  Approximating 
Function  Generator  (RMAFG)  developed  in  Task  4.  Table  1-3  lists  some  of  the 
potential  parameters  (independent  variables)  that  were  considered  in  the 
sensitivity  studies,  and  Table  1-4  lists  the  number  of  parameters  that  are 
associated  with  various  components  of  a typical  model  using  the  CVS  program. 

Results  from  Task  7,  Numerical  Verification  of  the  RMAFG,  are  also 
described  in  Section  4.  The  RMAFG  is  a multiparameter  polynomial  interpolating 
routine  which  is  useful  for  parametric  studies.  The  Fortran  listing  of  the 
RMAFG  program  is  presented  in  Appendix  B. 
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Table  l-l 


QUESTIONS  CONCERNING  REQUIRED  LEVEL  OF  MPPEL  PETAIL 

1.  Is  similitude  between  rigid  body  characterizations  of  selected  segments 
adequate  or  must  similitude  of  one  or  more  vibration  modes  also  be  preserved? 

2.  Are  simple  ball-and-socket  joints  an  adequate  characterization  of  the 
shoulder  or  must  the  separate  motions  of  the  clavicle  and  scapula  be  approxi- 
mated? 

3.  Is  it  sufficient  to  preserve  kinematic  similitudes  that  are  describable 
by  geometric  constraints  between  the  head  and  the  base  of  the  neck  or  must 
physical  or  computer  models  of  the  neck  include  dynamic  degrees  of  freedom  that 
correspond  to  the  human  neck? 

4.  Is  a planar  model  adequate  or  must  a three-dimensional  model  be  used  for 
conditions  essentially  symmetric  to  the  body? 

5.  Is  it  sufficient  to  preserve  similitude  of  impulses  and  coefficients 
of  restitution  during  "hard"  impact  processes  or  must  similitude  of  forces 
also  be  preserved? 

i 

6.  Is  it  necessary  to  account  for  the  non-symmetric  shape  of  joint  stops? 

7.  Is  a sliding/ rolling  characterization  of  a particular  contact  adequate 
or  must  similitude  of  the  deflection  characteristics  be  preserved? 

8.  Is  it  necessary  to  provide  separate  rigid  body  degrees  of  freedom  for 
the  hands  or  the  feet? 

9.  Must  the  similitude  of  surface  compliances  be  preserved  in  the  region 
of  belt  or  air  bag  contacts? 

10.  It  is  necessary  to  preserve  similitude  of  certain  dynamic  degrees  of 
freedom  which  include  inertial  effects  for  deformations  associated  with  contact? 

11.  With  respect  to  each  of  the  modeling  features,  and  for  various  typical 
crash  conditions  either: 

11.1  what  is  the  simplest  level  that  adequately  handles 
all  of  the  cases,  or 

11.2  what  is  the  most  detailed  model  that  is  within  the 
capabilities  of  the  CVS  program  and  that  is  inadequate 
for  one  or  more  of  the  cases? 

12.  What  is  the  simplest  and  least  expensive  mechanical  idealization, 
using  the  CVS,  that  appears  to  be  adequate  for  each  of  several  typical  crash 
conditions? 
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Table  1-2 


LEVELS  OF  SOPHISTICATION 


1.  Planar,  i.e.,  two-dimensional,  motion  using  planar  kinematic  constraint 

option  (for  symmetric  and  nearly  symmetric  crash  situations) . Include  a case 
where  the  two  sets  of  planes  formed  by  the  upper  and  lower  parts  of  the  arms 

and  legs,  respectively,  make  significant  angles  with  the  mid-sagittal  plane. 

2.  Abdomen  modeled  as  several  rigid  body  segments  each  with  independent  ro- 
tational degrees  of  freedom.  The  lumped  mass,  stiffness,  and  viscous  parameters 
should  be  consistent  with  the  Baseline  Model  (if  this  is  not  the  Baseline  Model 
itself)  . 

3.  Abdomen  modeled  the  same  as  for  item  (2)  except  that  the  individual 

abdomen  subsegments  have  orientations  with  respect  to  one  end  that  are  constrained 
to  be  functions  of  the  relative  orientations  of  the  opposite  end  of  the  abdomen. 

4.  Neck  modeled  similar  to  the  treatment  described  in  Item  (2)  for  the  abdomen. 

5.  Neck  modeled  similar  to  the  treatment  described  in  Item  (3)  for  the  abdomen. 

6.  Shoulder  model  that  includes  clavicle,  scapula  and  muscles  modeled  as  ten- 
sion elements  with  mass. 

7.  Globalgraphic  descriptions  of  joint  stops  that  correspond  to  the  human 
subjects . 

8.  Impulsive  contacts  for  nominal  conditions  involving  relatively  short 
duration,  i.e,  hard  contacts. 

9.  Degrees  of  freedom  that  approximate  certain  segment  modes  of  vibration 
during  events  that  include  hard  impacts.  (This  may  well  involve  modeling  seg- 
ments as  an  assembly  of  two  or  more  rigid  subsegments) . 
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Table  1-3 


LIST  OF  POTENTIAL  INDEPENDENT  VARIABLES  FOR  SENSITIVITY  STUDIES 


1. 

Components  of  inertia  tensors 

for  body  segments. 

2. 

Joint  friction  torque  unlock 

threshold. 

3. 

Joint  moving  friction  torque. 

4. 

Joint  stop  parameters. 

5. 

Flexural  moment  stiffnesses. 

6. 

Joint  locations. 

7. 

Initial  position  and  orientation. 

8. 

Contact  model  parameters. 
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Table  1-4 

SIGNIFICANT  PARAMETERS 


Significant  Parameters  No.  of  Parameters  Required/Element 


Segment 

mass  1 


inertia 

6 

contact  ellipsoid 

9 

integrator  tests 

1_2 

28 

TOTAL 

Ball  Joint 

location  and  orientation 

12 

spring  characteristics 

10 

viscous  characteristics 

7 

29 

TOTAL 

Pin  Joint 

location  and  orientation 

12 

spring  characteristics 

5 

viscous  characteristics 

7 

24 

TOTAL 

Euler  Joint 

location  and  orientation 

12 

spring  characteristics 

15 

viscous  characteristics 

11 

48 

TOTAL 

(Globalgraphic  joint  option  requires  a minimum  of  22  additional 
parameters) 


Belt 

reference  points 
slack 

stress-strain  zero  friction 
or  infinite  friction 


9 

1 

11  (mimimum) 
(22)  (minimum) 
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21  (32)  minimum  TOTAL 


Table  1-4  (continued) 


Airbag 

basic  parameters 

each  additional  reaction  panel 

Flexible  Element 

for  each  segment,  h function 

Contact  Surface 

size,  position  and  orientation 
force-deflection  characteristics 


38 

9 

38  + 9 x number  of  panels  TOTAL 
30 

30  x number  of  segments  TOTAL 
9 

50  typical  (minimum  16) 

59  TOTAL 
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Task  6,  the  Reduction  in  the  Number  of  Inertially  Related  Parameters, 
is  addressed  in  Section  5 where  a theorem  is  presented  which  shows  that  the 
mass  distribution  of  a set  of  connected  rigid  bodies  is  not  unique.  A sample 
case  is  presented  in  that  section,  and  the  Fortran  listing  of  the  computer 
program  for  computing  dynamically  equivalent  systems  is  contained  in 
Appendix  C. 
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CHECK  OF  THE  CVS  PROGRAM 


One  primary  purpose  of  Task  2 was  to  make  a complete  check  of  the  CVS 
program  and  the  algorithms  used  to  develop  the  program.  A problem  that  came  to 
light  was  the  failure  of  the  joint  constraints,  in  that,  in  some  runs  a noticeable 
drift  occurs  in  the  Pin  and  Euler  joints.  This  problem  was  studied  in  detail 
and  is  discussed  below. 

The  CVS  program  links  the  segments  by  imposing  forces  and  torques  of 
constraint  at  the  joints.  The  program  integrates  the  linear  accelerations 
and  velocities  of  the  reference  segment  (s)  to  obtain  the  linear  velocity  and 
the  position  of  these  segments.  It  then  uses  a chaining  algorithm  to  compute 
the  positions  and  the  velocities  of  the  other  segments.  The  use  of  this  chaining 
procedure  insures  consistency  of  the  linear  variables,  i.e.,  the  joints  do 
not  pull  apart.  However,  when  a joint  is  locked  or  pinned  (one  free  axis)  or 
has  one  locked  axis,  nothing  in  the  program  checks  or  corrects  for  any  inconsis- 
tencies that  may  occur  in  the  relative  angular  positions  or  velocities  because 
of  integration  or  precision  errors. 

To  state  the  problem:  if  segments  1 and  2 are  connected  by  a pin 

joint  (IPIN  = 1 or  IEULER  = 4,  5,  6),  the  angular  constraint  is  that  the  vector 
defining  the  pin  axis  in  segment  1 should  remain  parallel  to  the  vector  defining 
the  pin  axis  in  segment  2.  This  constraint  is  imposed  on  the  system  by  defining 
a torque  of  constraint  for  th.s  joint.  The  constraint  equation  is  obtained  in 
acceleration  form  by  twice  differentiating  the  equation  expressing  the  equality 
of  the  pin  axes.  Since  the  orientation  of  segment  1 is  updated  by  the  integrator 
independently  of  the  orientation  of  segment  2,  errors  in  integration  and  lack  of 
infinite  precision  in  the  calculations  will  cause  errors  in  orientation  (direction 
cosine  matrices)  and  hence  the  pin  vectors  may  not  remain  parallel.  This  drift 
effect  has  been  noticed  in  some  of  the  computer  runs.  For  simple  pin  joints  the 
maximum  error  detected  has  been  less  then  1 degree;  however,  for  the  Euler  joint 
with  one  locked  axis  a very  noticeable  drift  (several  degrees)  has  been  reported 
by  Leyland  Motors. 
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Since  the  typical  Part  572  dummy  joint  is  best  modeled  by  the  use  of 
the  Euler  joint,  any  drift  may  cause  significant  errors  in  the  simulation.  Thus, 
considerable  effort  was  expended  to  determine  the  cause  of  the  drift.  All  of 
the  subroutines  involving  the  joint  computations  were  examined  in  detail  to 
ensure  that  there  were  no  programming  errors  or  errors  in  the  algorithms.  No 
errors  were  found;  however,  it  should  be  noted  that  a correction  had  been  made 
to  subrouting  VISPR  in  1979.  This  correction  would  make  the  current  version  of 
the  program  perform  differently  from  previous  versions.  The  correction  did  not 
help  the  drift  problem. 

A simple  three  segment-two  joint  model  was  used  to  make  tests  on  the 
drift.  These  tests  verified  that  the  drift  was  a function  of  integrator 
accuracy.  As  the  ratio  test  on  the  integration  angular  acceleration  was  reduced, 
the  drift  was  reduced.  In  these  tests,  both  joints  were  Euler  joints;  one  had 
the  precession  axis  locked  and  the  other  had  the  spin  axis  locked.  It  was 
noticed  that  although  the  individual  angles  drifted  from  their  expected  values, 
the  sum  of  the  angles  was  close  to  the  expected  values.  The  reason  for  this 
has  not  been  determined. 

Modifications  to  subroutine  PDAUX  were  tried  in  an  attempt  to  correct 
the  drift.  In  particular  several  methods  of  correcting  the  quarternion  to  insure 
that  the  vector  part  of  the  quarternion  had  no  component  on  the  effective  locked 
axis  were  tried.  None  of  these  methods  has  a significant  effect  on  the  drift. 

It  appears  that  it  will  be  necessary  to  develop  a chaining  algorithm 
for  the  angular  positions  and  velocities  to  eliminate  the  drift  problem.  Until 
this  is  done,  it  is  recommended  that  the  integrator  tests  be  used  to  hold  the 
drift  to  an  acceptable  level.  No  method  of  predetermining  the  best  values  to 
use  for  these  tests  is  known;  hence  they  will  have  to  be  determined  by  trial 
and  error.  The  tabular  time  histories  of  the  joint  behavior  show  the  values 
of  the  angles  involved.  Examination  of  these  time  histories  allows  the  user  to 
detect  the  drift  in  any  angle  that  should  be  constant  (locked  axis),  but 
unfortunately  the  drift  or  error  in  a time-varying  angle  cannot  be  determined. 
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3.0 


COMPARATIVE;  STUDY  OF  LEVELS  OF  SOPHISTICATION 


A significant  feature  of  the  CVS  program  is  the  ability  to  vary  the 
complexity  of  the  model  of  a physical  system.  For  most  applications  simulating 
the  Part  572  dummy,  a fifteen- segment  model  has  been  used.  However,  the 
skeletal  structure  of  the  Part  572  dummy  consists  principally  of  33  rigid 
segments  plus  the  rubber  neck  and  the  rubber  lumbar  spine.  Vinyl  foam  "flesh" 
surrounds  some  of  these  segments  to  form  an  approximate  human  shape.  No  internal 
organs  are  included  except  for  a visceral  sac  in  the  lower  torso. 

Since  the  CVS  program  requires  that  the  structure  be  defined  as  a 
set  of  rigid  segments  connected  by  joints,  the  neck  and  the  spine  each  must  be 
approximated  as  a set  of  rigid  segments  and  the  "flesh"  is  assumed  to  be  rigid 
and  part  of  the  segments  to  which  it  is  attached.  This  study  did  not  address 
the  effects  of  this  assumption  for  the  "flesh."  The  spring  and  viscous 
characteristics  attributed  to  the  joints  are  primarily  due  to  the  interference 
between  the  flesh  of  the  segments  connected  by  each  joint.  In  studies  performed 
in  another  program  (Reference  1)  it  was  found  that  the  visceral  sac  significantly 
affected  the  dynamic  motion  of  the  torso  and  had  to  be  accounted  for  in  defining 
the  properties  of  the  joints  representing  the  lumbar  spine. 

Many  of  the  33  rigid  segments  of  the  dummy  are  components  of  the  joints 
which  are  modeled  by  using  the  "Euler"  joint  routine  in  the  program.  This  routine 
accurately  models  the  kinematics  of  the  joint  but  ignores  the  inertial  reactions. 
The  inertial  reactions  are  approximated  by  including  the  mass  and  inertia  of 
the  joint  components  as  part  of  the  adjoining  segments.  For  example,  the  elbows 
are  modeled  as  "Euler"  joints  and  the  mass  of  the  joint  structure  is  included 
in  the  masses  of  the  upper  and  lower  arms.  Thus,  by  use  of  the  "Euler"  joint 
routine  for  the  elbows,  ankles,  and  hip  joints,  six  of  the  segments  were  combined 
with  others  to  reduce  the  complexity  of  the  model  to  27  segments.  The  knees 
were  modeled  as  pin  joints  and  the  upper  legs  as  single  segments.  In  addition, 
the  hands  were  included  as  part  of  the  lower  arms,  further  reducing  the  complexity 
to  21  segments.  (In  some  applications  it  may  be  desirable  to  allow  for  motion 
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of  the  wrist -j  Each  shoulder  structure  of  the  Part  572  dummy  includes  three 
segments.  Replacement  of  these  by  a single  "Euler"  joint  for  each  shoulder 
reduces  the  complexity  to  15  segments.  A significant  part  of  the  effort  on 
this  contract  was  directed  toward  a study  of  the  shoulder  which  is  discussed 
in  a following  section. 

3 . 1 General  Discussion  of  the  Levels  of  Sophistication 

* i 

Various  degrees  of  complexity  in  modeling  certain  of  the  dummy 
characteristics,  in  particular  the  rubber  neck,  the  rubber  lumbar  spine  (abdomen) 
and  the  shoulder  were  investigated.  These  degrees  of  complexity  are  referred 
to  as  "levels  of  sophistication"  and  are  summarized  in  Table  1-2.  The  main 
questions  regarding  the  required  level  of  model  detail  are  given  in  Table  1-1. 

Level  of  sophistication  1 of  Table  1-2,  the  study  of  planar  motion, 
refers  to  the  symmetry  options  available  in  the  CVS  program.  There  are  two 
symmetry  options  in  the  program.  The  first  is  the  two-dimensional  symmetry 
option  whereby  all  motion  perpendicular  to  the  mid-sagittal  plane  is  suppressed. 
The  second  option  provides  mirror  symmetry  where  the  motion  is  reflected  in  the 
mid-sagittal  plane,  e.g.,  if  the  right  arms  and  legs  move  to  the  right,  the  left 
arms  and  legs  will  move  the  same  amount  to  the  left.  The  computer  runs  that  were 
made  using  these  options  showed,  in  some  cases,  a reduction  in  computer  time  when 
compared  to  the  same  cases  run  without  the  symmetry  options.  This  is  primarily 
due  to  the  smoothness  of  the  integration  which  in  turn  allowed  the  use  of 
larger  integration  steps  and,  hence,  fewer  computations.  The  user  is  advised  to 
use  these  options  whenever  they  are  applicable  to  his  problem. 

.1 

Levels  of  sophistication  2,  3,  4,  and  5 of  Table  1-2  are  concerned 
with  the  modeling  of  the  dummy  spine  and  neck.  The  15  segment  Baseline  model  of 
the  Part  572  dummy  (Ref.  1,  Volume  2)  uses  a single  rigid  segment  to  model  the 
abdomen  (lumbar-spine)  and  a single  rigid  segment  to  model  the  neck.  Both  of 
these  segments  are  connected  to  the  adjoining  segments  with  ball  joints  in  the 
model . 


The  Head/Neck  Pendulum  tests  that  are  reported  in  Reference  1, 

Volume  2 were  studied  in  detail.  It  was  apparent  from  these  results  that  at 
least  one  segment  must  be  used  for  the  neck.  For  example,  the  behaviour  of 
the  angle  of  the  head  with  respect  to  the  pendulum  in  the  first  ten 
milliseconds  as  shown  in  Figure  4-3  of  the  cited  reference  could  not  be 
obtained  with  a model  in  which  the  head  was  connected  directly  to  the  upper 
torso.  Results  of  simulations  in  which  the  neck  was  assumed  to  be  composed 
of  either  two  or  three  segments  showed  no  significant  differences  in  response 
from  the  one  segment  model.  The  Head/Neck  Pendulum  tests  also  indicated  some 
extension  (stretching)  of  the  neck.  Attempts  to  simulate  this  effect  by 
replacing  the  ball  joints  with  sets  of  springs,  thus  effectively  disconnecting 
the  neck  from  the  upper  torso  and  the  head,  were  unsuccessful.  The  failure 
was  due  primarily  to  the  inability,  using  the  current  capabilities  of  the  CVS 
program,  to  define  a set  of  springs  that  would  allow  extension  but  would  not 
allow  significant  lateral  motion  at  the  joint.  This  deficiency  could  be 
removed  by  adding  a slip-joint  capability  to  the  CVS  program.  Such  a joint  could 
be  defined  as  one  imposing  a kinematic  constraint  that  would  allow  relative 
linear  motion  at  the  joint  parallel  to  fixed  lines  in  each  of  the  adjoining 
segments . 


Simulations  of  torso  pendulum  impact  tests  reported  in  Reference  1, 
Volume  2 indicated  that  the  principal  deficiency  of  the  single-segment  model 
of  the  abdomen  was  not  in  the  modeling  of  the  spine  but  in  the  lack  of  modeling 
of  the  visceral  sac  that  is  in  the  Part  572  dummy.  No  significant  stretching  of 
the  dummy  abdomen  occurs  because  of  the  constraint  of  the  steel  cable  that 
passes  through  the  center  of  ijhe  rubber  spine.  More  detailed  models  of  the 
abdomen  were  not  studied  in  the  current  effort. 

Level  of  sophistication  6,  as  stated,  applies  to  the  human  victim 
only  since  the  shoulder  of  the  dummy  is  a mechanical  linkage.  The  shoulder 
linkage  of  the  Part  572  dummy  is  studied  in  some  detail  in  the  next  section 
of  this  report.  The  consideration  of  a model  for  the  human  shoulder  was  not 
within  the  scope  of  this  research  program. 
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Level  of  sophistication  7,  dealing  with  the  g lobalgraphic  description 
of  joint  stops  is  applicable  whenever  there  is  sufficient  motion  for  a joint  to 
enter  a stop.  The  user  is  advised  to  specify  this  option  whenever  there  is  a 
probability  of  a joint  entering  a stop. 

Level  of  sophistication  8,  impulsive  contacts,  can  be  examined  with 
a simpler  model  than  that  of  the  entire  dummy.  Question  number  5 of  Table  1-1 
is:  Under  what  conditions  is  it  possible  to  replace  a contact  with  a hard 

surface  with  an  idealized  impulse?  The  use  of  an  impulse  is  justified  when 
there  is  a change  in  velocity  in  a short  enough  time  interval  that  there  is 
no  significant  change  in  position.  The  problem  is  the  definition  of  a 
significant  change  in  position.  Section  6.2  of  this  report  presents  a more 
detailed  discussion  of  this  problem. 

The  use  of  the  impulse  in  relation  to  question  10  of  Table  1-1 
should  not  be  overlooked.  The  program  has  an  inertial  spike  option  which 
models  some  of  the  inertia  effects  of  a contact.  It  was  originally  developed 
for  the  case  of  a head  striking  a windshield  where  the  mass  of  glass  must  be 
accelerated  when  it  shatters.  This  inertial  effect  usually  lasts  only  a few 
milliseconds.  The  transfer  of  energy  could  be  approximated  by  use  of  an 
impulse  with  a negative  coefficient  of  restitution  (the  program  allows  this) 
which  would  allow  the  impacting  object  to  continue  in  the  same  line  of  motion 
without  a reversal  of  velocity. 

Level  of  sophistication  9 is  concerned  with  transverse  modes  of 
vibration  which  may  be  inducec,  when  a hard  contact  is  made  with  a segment 
such  as  a leg  and  the  leg  bends  due  to  the  severity  of  the  impact.  This  may 
be  modeled  by  use  of  several  subsegments  but  the  stiffness  of  the  joints  may 
require  a small  integrating  step  and  hence  make  it  economically  unfeasible. 

The  flexible  element  option  may  alleviate  this  difficulty.  The  long  bones  in 
the  dummy  are  steel  shafts  and  hence  are  much  stiffer  than  a human  bone. 

Thus,  segment  transverse  vibrations  may  be  a more  important  consideration  for 
humans  than  for  a dummy  model,  The  consideration  of  longitudinal  modes  of 
vibration  are  beyond  the  capabilities  of  the  CVS  program. 
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3.2 


Study  of  the  Part  572  Dummy  Shoulder 


In  the  Baseline  model  of  the  Part  572  dummy  (Reference  1),  the 
shoulder  is  represented  as  an  Euler  joint  with  the  spin  axis  locked.  The 
actual  mechanical  structure  is  illustrated  in  Figure  3-1  and  consists  of  three 
principal  pieces  for  each  shoulder.  These  are  the  stemo-clavicular  link,  the 
clavicle,  and  the  shoulder  yoke.  The  stemo-clavicular  link  pivots  about  an 
axis  which  is  tilted  (pitched)  5 degrees  to  the  x (forward)  axis  of  the  upper 
torso.  The  motion  is  constrained  by  a highly  nonlinear  viscous  damping 
mechanism  and  by  hard  mechanical  stops  which  limit  the  elevation-depression 
motion  of  the  shoulder  girdle  to  15  degrees.  The  clavicle  is  pinned  to 
this  link  with  the  pin  parallel  to  the  nominal  z (down)  axis  of  the  torso. 

The  motion  about  this  pin  is  also  restricted  to  approximately  plus  and  minus 
15  degrees.  The  shoulder  yoke  connects  the  upper  arm  to  the  clavicle.  At 
the  clavicle  end  is  a pin  joii^t  approximately  parallel  to  the  y axis  of  the 
torso,  and  at  the  upper  arm  end  of  the  yoke  is  another  pin  joint  whose 
axis  is  approximately  parallel  to  the  x axis  of  the  torso.  These  joints  allow 
flexion-extension  and  abduct ion- adduct ion  motion  of  the  upper  arm,  respectively. 

In  the  Baseline  model,  the  shoulder  yoke  is  modeled  by  an  Euler 
joint  and  the  sterno-clavicular  link  and  the  clavicle  are  not  modeled.  The 
purpose  of  this  study  was  to  examine  the  effects  of  using  a shoulder  model 
that  more  closely  represents  the  actual  dummy  shoulder  assembly. 

The  computer  model  for  investigating  the  behavior  of  the  shoulder 
was  based  on  the  model  used  to  simulate  the  torso-pendulum  tests  described 
in  Reference  1.  In  the  torso-pendulum  tests  the  pendulum  is  raised  to  an 
almost  horizontal  position  and,  when  released,  strikes  a honeycomb  structure 
which  stops  the  motion  of  the  pendulum  when  it  reaches  a vertical  orientation. 

The  torso  and  associated  segments  are  thus  subjected  to  an  impulsive  force  at 
the  point  where  the  lower  torso  is  attached  to  the  pendulum.  The  three-piece 
shoulder  structures  and  the  upper  and  lower  arms  were  added  to  the  torso-pendulum 
test  simulation  model  to  form  a fourteen-segment,  fourteen- j oint  model.  The 
segments  and  joints  are  listed  in  Table  3-1  and  the  configuration  is  illustrated 
in  Figure  3-2. 
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SEGMENTS  AND  JOINTS  OF  MODEL  USED  FOR  STUDY  OF  SHOULDER 
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Figure  3-2  MODEL  OF  PENDULUM  TEST  OF  TORSO  WITH  SHOULDER  AND  ARMS 
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A series  of  computer  simulation  runs  was  made  in  which  selected 
sets  of  the  shoulder  joints  were  locked  or  unlocked  for  three  different 
initial  pendulum  impact  velocities  and  three  different  orientations  of  the 
torso-shoulder-arm  model  on  the  pendulum.  A particular  set  of  six  runs  that 
was  made  is  tabulated  below: 


Joint 

Status 

Level  of  Sophistication 

Run 

X 

Z 

Y 

P 

N 

s 

1 

U 

u 

U 

L 

U 

L 

Most  sophisticated  - all  segments  active 

2 

U 

u 

L 

U 

U 

L 

Sterno-clavicle  free-Euler  for  shoulder 

3 

U 

L 

L 

U 

u 

L 

Clavicle  locked-Euler  for  shoulder 

4 

L 

U 

L 

U 

u 

L 

Clavicle  free-Euler  for  shoulder 

5 

L 

L 

L 

U 

u 

L 

Least  sophisticated-CVS  baseline 

6 

L 

L 

U 

L 

u 

L 

Sterno-clavicle  locked,  shoulder  yoke- 
clavicle  pivot  used  instead  of  Euler 

where  L indicates  that  the  joint  is  always  locked,  and 

U indicates  that  the  joint  is  allowed  to  unlock 
at  a specified  torque. 

The  column  labels  for  joint  status  are  (see  Table  3-1  and  Figure  3-1): 

X for  RX  and  LX,  the  joint  connecting  the  stemo-clavicle  link 
to  the  upper  torso, 

Z for  RZ  and  LZ,  th  i joint  connecting  the  stemo-clavicle  link 
to  the  clavicle, 

Y for  RY  and  LY,  th<?  joint  connecting  the  clavicle  to  the 
shoulder  yoke. 

P the  precession  axis  of  the  Euler  joint.  This  joint  is  a 
direct  replacement  for  Y;  when  one  is  locked  the  other 
should  be  unlocked, 

N the  nutation  axis  of  the  Euler  joint  which  should  never  be 
permanently  locked. 

S the  spin  axis  of  the  Euler  joint.  This  joint  is  always 
permanently  locked. 
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The  initial  angular  velocity  of  the  pendulum  and  all  of  the  segments  was  113 
degrees  per  second  about  the  y axis  (pitch).  This  represents  an  initial 
linear  velocity  of  148  inches  per  second  at  the  eg  of  the  lower  torso  and  184 
inches  per  second  at  the  shoulder.  The  arms  were  initially  parallel  to  the  z 
axis  of  the  torso.  All  joints  except  2,  3 and  14  were  initially  locked. 

Plots  of  selected  response  measures  for  each  of  the  six  cases  are 
contained  in  Appendix  A.  The  response  measures  plotted  are: 

(a)  The  resultant  linear  acceleration  of  the  right  lower  arm  versus  time. 

(b)  The  z versus  x displacement  of  the  right  lower  arm. 

(c)  The  y versus  x displacement  of  the  right  lower  arm. 

(d)  The  flexure  angle  of  the  shoulder  yoke-clavicle  pivot  joint  versus  time. 

(e)  The  precession  and  nutation  of  the  Euler  joint  at  the  shoulder  versus  time. 

S ■ 

(f)  The  precession  and  nutation  of  the  Euler  joint  at  the  elbow  versus  time. 

(g)  The  resultant  linear  acceleration  of  the  upper  torso  versus  time. 

(h)  The  pitch  of  the  upper  torso  versus  time. 

The  results  obtained  with  the  most  sophisticated  model  (Run  1)  and 
with  the  least  sophisticated  model  (Run  5)  of  the  shoulder  are  overlaid  in 

1 

Figure  3-3  for  comparison.  It  may  be  noted  in  Figure  3-3 (d)  that  the  flexure 
angle  at  the  shoulder  when  the  shoulder-clavicle  is  free  (Run  1 with  Y unlocked) 
seems  to  reverse  direction  at  approximately  85  milliseconds.  This  is  because 
the  shoulder-clavicle  pivot  is  modeled  as  a pin  joint  which  is  not  direction 
sensitive.  The  first  part  of  the  curve  should  be  reflected  about  the  abscissa 
so  the  plot  would  be  much  like  that  of  the  precession  angle  of  the  shoulder 
for  Run  5 (precession  axis  of  the  Euler  joint  unlocked)  shown  in  Figure  3-3(e). 

No  definitive  conclusion  was  reached  from  the  study  except  to  note  that 
the  response  of  the  upper  torso  is  essentially  the  same  in  all  of  the  runs.  The 
response  of  the  lower  arms  is  different  in  each  of  the  runs.  The  reader  should 
draw  his  own  conclusions  as  to  whether  or  not  the  differences  are  significant. 
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Figure  3-3  COMPARISON  OF  RESPONSES  FROM  TWO  SHOULDER  MODELS 
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Figure  3-3  (Cont'd. 
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(f)  RE  PRECESSION  AND  NUTRTIGN  VS.  TIME 
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(g)  UPPER  TORSO  RESULTANT  ACCELERATION  VS.  TIME 
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ch)  UPPER  TORSO  PITCH  VS.  TIME 


4.0 


RESPONSE  MEASURE  APPROXIMATING  FUNCTION  GENERATOR 


In  some  applications  of  the  CVS  program  the  user  may  wish  to 
perform  various  parameter  studies  in  which  he  varies  one  or  more  parameters 
in  the  model  and  notes  the  changes  in  selected  response  measures.  In 
practice,  the  resulting  response  measures  are  usually  plotted  and  smooth 
curves  drawn  to  show  the  variation  with  the  parameters.  It  is  natural  to 
assume  that  the  value  of  the  response  measure  for  any  intermediate  value  of 
a parameter  can  be  obtained  from  the  curve,  i.e.,  the  curve  is  used  as  an 
interpolating  function. 

The  use  of  a hand  drawn  curve  is  quite  practical  if  only  one 
parameter  has  been  varied,  but  if  two  parameters  are  varied  a surface  must 
be  drawn  to  use  this  graphical  interpolation  approach.  If  more  than  two 
parameters  are  varied,  the  interpolation  problem  usually  becomes  too 
complicated  to  solve  using  planar  graphs. 

Hence,  a multiparameter  interpolating  routine  was  developed  to 
assist  the  users.  This  routine  is  called  the  Response  Measure  Approximating 
Function  Generator  (RMAFG).  The  mathematical  development  and  a description 
of  the  Fortran  program  of  the  RMAFG  is  contained  in  the  following  sections 
along  with  a sample  application. 


29 


4.1 


Mathematical  Development 


Let  y be  a scalar  function  of  the  N parameters  X.  3 J = 13  n3 


i.e.,  y = f (Xj3  Xg  ...3  X^) . Assume  that  this  function  may  be  approximated 
as  a polynomial  in  the  N parameters  as: 


y - c 


o 


N 

+ Z 
i=l 


C.  X. 


N 

+ Z 


N 

l 


i=l  j=l 


C..  X.  X. 

T'J  1 J 


N 

+ l 


N 

Z 


N 

Z 


i=l  j=l  k=l 


+ etc . 


C . X.  X.  xn 

k ^ j A: 


Constant  Linear 


Linear 


Quadratic 


Cubic 


Quartic 3 etc. 


(4.1) 


The  number  of  terms  of  each  degree  is  given  by  the  binomial  coefficient 


I 


,N-1+K , 
K 


C t™)  where  K is  the  degree.  The  total  number  of  terms  up  to  and  including 


,N+K 


degree  K is  given  by  f ^ ) . This  is  tabulated  in  Table  4-1  below. 


Table  4-1 


NUMBER  OF  TERMS  IN  MULTINOMIAL  EXPANSION 


! 


No.  of  Parameters  Degree  K 


N 

2_ 

3 

4 

5 

6 

7_ 

1 

2 

3 

4 

5 

6 

7 

8 

2 

3 

6 

10 

15 

21 

28 

36 

3 

4 

10 

20 

35 

56 

84 

120 

4 

5 

15 

35 

70 

126 

210 

330 

5 

6 

21 

56 

126 

252 

462 

792 

6 

7 

28 

84 

210 

462 

924 

1716 

7 

8 

36 

120 

330 

792 

1716 

3432 
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Note  that  the  first  row  and  the  first  column  of  the  table  are  just 
a sequence  of  consecutive  integers  and  that  the  other  entries  are  the  sum 

of  the  entry  above  and  the  entry  to  the  left,  i.e.,  a.  .=  a.  7 . + a.  . 7. 

T' 3 J J 1 

In  a typical  sensitivity  study  some  set  of  parameters,  X.  3 is 

J 

varied  over  a preselected  set  of  values  and  the  values  of  a response  measure, 
y j are  recorded.  The  values  of  the  response  measure  may  then  be  plotted  as 
functions  of  the  various  variables.  However,  if  one  is  interested  in 
estimating  the  value  of  the  response  measure  for  a set  of  parameters 
different  from  one  of  the  preselected  sets  it  is  necessary  to  either  compute 
y for  the  new  set  (by  running  the  program)  or  to  estimate  y by  inter- 
polating  the  data  available.  A functional  representation  such  as  given  by 
Equation  4.1  provides  a means  of  interpolation.  The  coefficients,  C , can  be 
calculated  from  the  values  of  y obtained  from  the  preselected  sets  of  X_.  rs 
by  a least  square  procedure.  The  number  of  independent  sets  of  data  must  be 
at  least  equal  to  the  number  of  coefficients  as  given  in  Table  4-1  but  a 
significantly  greater  number  than  this  minimum  is  recommended  to  achieve  a 
certain  degree  of  smoothing  from  the  least  square  procedure. 

It  is  difficult,  if  not  impossible,  to  cite  how  many  values  should 
be  used  because  of  the  complexity  of  the  problem.  It  is  easier  to  draw  con- 
clusions 'after  the  fact';  that  is,  the  coefficients  are  computed  for  a given 
set  of  data  and  the  function  evaluated  for  intermediate  points.  If  there  are 
rapid  variations  of  the  function,  its  use  as  an  interpolation  function  is  in 
doubt  whereas  if  the  function  behaves  smoothly,  then  it  may  be  useful.  For 
example,  consider  a cubic  fit  to  four  points  of  a function  of  one  variable. 
Both  a rapid  variation  and  a smooth  variation  are  illustrated  in  Figure  4-1. 

In  both  cases  the  cubic  goes  through  all  of  the  data  points  but  the  rapidly 
varying  function  clearly  can  result  in  less  accurate  interpolated  values. 


Poor  Fit  - Rapidly  Varying 
Function 


Good  Fit  - Smoothly  Varying 
Function 


Figure  4-1  EXAMPLES  OF  POOR  AND  GOOD  FITS  TO  DATA  USING  CUBIC  FUNCTIONS 


Returning  to  Equation  4.1,  if  the  value  of  each  X is  replaced  by 

its  deviation  from  some  nominal  value,  i.e.,  replace  X.  by  X.  - X . 3 X. 

Ts  'is  'Ll)  J 

by  Xj  - % -Q,  etc.,  then  the  coefficients  in  the  expansions  may  be  considered 
as  the  terms  in  a Taylor  series.  That  is,  C is  the  value  of  the  function  at 
the  nominal  point  X , i = 1,N  and  the  linear  terms  C\  are  the  partial 
derivatives  evaluated  at  the  nominal  point,  etc. 


4.2 


RMAFG  Fortran  Program 


A Fortran  program  has  been  developed  to  evaluate  the  coefficients  in 
the  expansion  of  the  form  givun  by  Equation  4.1.  It  consists  of  a Main 
Program  and  four  subroutines.  The  basic  coding  is  written  in  a fashion  such 
that  the  size  of  the  problem  Is  limited  only  by  the  storage  assignments  made 
in  the  dimension  statements.  In  the  listing  of  the  program  given  in  Appendix  B, 
the  program  is  dimensioned  to  handle  ten  independent  variables  (X.)  and  100 

'Is 

data  points. 

Three  modes  of  operation  are  provided: 
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! I 


Mode  0 --  the  expression  as  given  by  Equation  4.1  is  treated. 


Mode  1 --  the  X.  of  Equation  4.1  are  replaced  by  X. 

'Is  ^ 

is  the  average  of  the  given  values,  i.e.. 


X.n  where  X. 

1*0  zO 


Xi0  M 


M 

l 

k=l 


X.(k )a  i = 13N 
z 


X.(k)  is  the  value  of  X.  in  the  kth  data  set. 

z z 

Mode  2 --  The  X^  of  Equation  4.1  are  replaced  by  (X^  - X.^)/o 

where  X.n  is  as  in  Mode  1 and 

zO 

M 

1 2 21/2 

a = (—  2 X.(k)  - X.J  .is  the  standard  deviation. 

M i z zO 

Mode  2 has  an  advantage  of  normalizing  the  parameters  (independent 
variables)  to  a common  distribution  so  that  the  magnitude  of  the  C’s  are 
more  readily  comparable. 

The  program  allows  the  degree  of  parameter  to  be  limited.  For 

example,  we  may  have  3 parameters  X^3  X ^ and  X ^ and  use  a cubit  fit  on  all 

parameters,  in  which  case  all  20  coefficients  ( C’s ) will  be  evaluated;  or, 

we  may  limit  X j to  the  1st  degree,  x % to  the  2nd  degree  and  X ^ to  the 

3rd  degree.  In  this  case  the  terms  will  be  y - Cr  + C-iX-i  + Wn 

o P ~U  1 1 Z 2 

+ + ?22  X2  + CZ2X2X"5  + CZ^Z  + CZZ^Z’  The  only  restriction  of 

the  program  is  that  the  parameters  be  labeled  so  that  they  are  ordered  on 

increasing  degrees,  as  in  the  example  just  cited. 

The  mathematical  procedure  is  most  readily  described  by  treating  the 
C’s  and  the  combinations  of  the  parameters  as  components  of  vectors.  In  the 
above  example  let: 


and 


r~  = (c  r r c r r r r > 

- ! O3  l3  ^23  Z3  22 3 323  333  333 


i?  — fly  y y y ^ 

- 3 V 23  3 3 23  “2“33  n 3 3 " 3 


x0x7J  x2  xh. 


Equation  4.1  may  be  written  as: 


y 


m 


If  we  use  the  subscript  m to  denote  the  mth  data  point  then 


The  error  of  the  fit,  E 3 is  then 


m 


E = y 
m 3 m 


- iFc 


- m — 


The  mean  squared  error  is 


* - i " 

m - m 


m 


M 


(y  - CTU  ) (u  - l/c) 

Jm m sm  -m— 


= y2  + (?  UUFC_  - C?  Uy  - ylF  C 


where 

T 
UU 


yv 


m 

Z U 

rf  , 

an 

m 

m=l 

-m 

A 

m 

m 

E y 

m 
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m 
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The  mean  squared  error  is  minimized  when  C has  the  value 
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C = (rjr/)  1 Uy 


squared  error  is: 


With  ,this  value  of  C_  the  average  error  E is  zero  and  the  average 


E2  = y2  - ylF  C. 

A special  matrix  inversion  routine  was  written  to  compute  the  C. 
It  allows  the  program  to  output  the  results  obtained  by  the  lower  order 
functions.  That  is,  if  a cubic  fit  is  specified  it  evaluates  the  best 
constant,  the  best  linear  fit,  the  best  quadratic  fit,  and  finally,  the 
best  cubic  fit. 


The  mean  squared  erfor  is  computed  for  each  of  these  intermediate 
fits  and  is  printed  along  with  the  values  of  the  coefficient  C_. 

The  general  flow  of  the  program  is  as  follows. 

Read  N3  (NDG(J),j=lsN) 

N number  of  parameters  (independent  variables) , 

if  N = 0 program  will  terminate 

NDG(J)  maximum  degree  of  parameter  X^. 

These  must  be  ordered:  NDG(J-l)  <_NDG(J). 

Call  SET  Z. 


Subroutine  SET  Z computes  the  subscripts  used  to 
identify  the  C's  and  checks  whether  any  storage 
limits  are  violated. 
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Read  NDATA,  MODE,  IVAL. 


NDATA,  number  of  data  sets. 

MODE,  0 normal 

1 remove  mean 

2 remove  mean,  normalize  by  a , (the  standard  deviation). 


IVAL,  0 

1 


compute  fit  and  go  back  for  another  case, 

! 

read  in  additional  data  and  evaluate  fit 
at  these  points. 


Read  'Y(J),  (X(J,I)  (J=1,N) , 1=1,  NDATA. 

\ 

Y(I)  function  value  at  point  I 


X(J,I)  parameter  value  at  point  I 
Call  SIGN. 


Subroutine  SIGM  evaluates  the  mean  and  sigma  of  y 
It  modifies  X.  depending  on  the  mode  selected. 


and  the  X .. 

J 


T — 

The  program  then  computes  the  UU  matrix  and  the  Uy  vector  and 
calls  subroutine  SETU  which  evaluates  the  U vector  for  each  data  point. 

T 

Next,  subroutine  SOLVE  is  called  which  inverts  the  UU  matrix  and 
solves  for  the  C's.  Subroutine  SOLVE  uses  a bordering  scheme  so  that  the 
intermediate  results  are  available.  As  each  fit  is  completed  (constant,  linear, 
quadratic,  etc.),  the  main  routine  evaluates  the  error  of  the  fit  and  outputs 
the  error  and  the  values  of  the  coefficients.  The  routine  checks  IVAL  to  see 
if  further  evaluations  are  desired.  If  so,  it  reads  the  desired  values  of  the 
parameter  and  evaluates  y. 


4.3 


Sample  Numerical  Evaluation  of  the  RMAFG 


Using  the  Head-Neck  Pendulum  test  simulation  (Reference  1, 

Volume  2)  as  a baseline  run,  the  linear  spring  coefficients  of  the  neck  joints 
were  varied  from  28  to  34.4  in  steps  of  1.6,  the  viscous  friction  of  the  neck 
joints  was  varied  from  0.130  to  0.170  in  steps  of  0.01  and  the  Coulomb 
friction  of  the  neck  joints  was  varied  from  90  to  110  in  steps  of  5.  This 
produced  a series  of  125  runs  (each  variable  has  five  different  values). 
Several  response  measures  were  extracted  from  the  tabular  time  histories 
and  processed  by  the  RMAFG.  The  basic  results  are  summarized  in  Table  4-2. 

Table  4-2 

SELECTED  OUTPUT  OF  THE  RMAFG 

Root  Mean  Square  Error  For 
Degree  of  Fit 


Response  Measure 

Mean 

Sigma 

Constant 

Linear 

Quadratic 

Cubic 

1.  HIC 

91. SI 

1.65 

1.6524 

0.4662 

0.2697 

0.2052 

2.  HSI 

122.85 

3.11 

3.1068 

0.4241 

0.3342 

0.2971 

3.  Time  of  zero  head  pitch,  ms. 

103.13 

3.16 

3.1570 

0.1522 

0.0135 

0.0044 

4.  Time  of  max.  head  resultant 
acceleration,  ms. 

33.01 

0.05 

0.0533 

0.0430 

0.0275 

0.0151 

5.  Value  of  max.  head 

resultant  acceleration,  g 

29.88 

0.24 

0.2352 

0.1551 

0. 1151 

0.0750 

6.  Time  of  max.  head  pitch, 
ms . 

55.94 

1.17 

1.1749 

0.0590 

0. 0221 

0.0160 

7.  Value  of  max.  head  pitch. 

61.95 

2.26 

0.2600 

0.1099 

0.0249 

0.0239 

deg. 


The  time  of  zero  head  pitch  was  estimated  by  linearly  interpolating 
the  values  of  pitch  just  before  and  just  after  the  pitch  passed  through  zero. 
In  all  runs,  the  tabular  time  histories  were  obtained  at  2 millisecond 
intervals.  The  times  of  maximums  and  the  values  of  the  maximums  were  obtained 
by  a quadratic  interpolation  formula  using  the  three  points  from  the  tabular 
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time  histories  in  the  vicinity  of  the  maximum.  This  quadratic  interpolation 
formula  is  -- 

f = fz  + lWf3  - fjXt  - t2)/d  + l/2(f:  - 2f2  + fg)(t  - t2)2/dZ 
where  f is  the  value  of  the  function  at  t and  d = 

The  time  of  the  maximum  is  given  by 
and  the  value  of  the  maximum  is 

In  using  the  RMAFG  the  three  parameters  were  identified  as  follows: 
Xj  linear  spring  coefficient 

Xg  viscous  friction 

x ^ coulomb  friction. 

The  mean,  m , and  the  standard  deviation,  a , of  each  of  the 
parameter  variations  were  computed  and  the  normalized  variables  (mode  2)  were 
used  by  the  RMAFG,  i.e.,  the  variables  used  were  Z.  = (x  . - m) /o . The  use 

^ Is 

of  these  normalized  variables  has  the  advantage  of  allowing  one  to  directly 
compare  the  significance  of  the  coefficients  of  the  fit  obtained  by  the  RMAFG. 

With  reference  to  Table  4-2,  which  shows  the  root  mean  square  error 
of  the  fit  as  a function  of  the  degree  of  the  fit,  we  note  the  following: 
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1.  For  all  responses  except  4 and  5 which  deal  with  the 
head  acceleration  there  is  a dramatic  reduction  in 
the  rms  error  when  the  degree  is  changed  from  a 
constant  to  a linear  fit. 

2.  The  most  significant  improvement  resulting  from  increasing 
the  degree  from  a linear  to  a quadratic  fit  occurs  for 
responses  3 and  7. 

3.  Response  3 shows  the  most  significant  improvement 
resulting  from  changing  from  a quadratic  to  a cubic  fit. 

One  conclusion  that  may  be  drawn  is  that,  for  response  7 (the  value 
of  maximum  head  pitch) , a quadratic  fit  is  sufficient  since  use  of  a cubic  . 
results  in  a negligible  reduction  of  the  error.  However,  note  that  the 
standard  deviation  of  the  maximum  pitch  angle  is  2.26  degrees  and  the  mean 
value  is  61.95  degrees  (computed  from  the  125  runs).  Thus,  even  though  a 
quadratic  fit  reduces  the  error  by  a factor  of  four  from  the  linear  fit,  the 
linear  fit  approximates  the  data  with  a rms  error  of  only  0.11  which  may  be 
adequate  for  most  purposes. 

From  a practical  point  of  view,  there  is  not  much  variation  of  any 
of  the  response  measures  selected  for  this  set  of  runs.  The  most  significant 
deviations  are  those  related  to  times,  responses  3 and  6,  and  the  value  of 
maximum  head  pitch,  response  7.  The  use  of  a linear  fit  in  these  three  cases 
reduces  the  error  to  a small  fraction  of  the  standard  deviation. 

It  is  recommended  that  the  user  use  the  RMAFG  to  evaluate  the  error 
reduction  that  results  from  increasing  the  degree  of  the  function  so  that  the 
degree  of  fit  that  is  adequate  for  his  purpose  may  be  selected.  It  is 
further  recommended  that  interpolated  results  at  intermediate  points  (points 
other  than  the  inputted  data  used  for  the  fit)  be  checked  to  ensure  that  the 
degree  selected  does  not  give  a 'poor'  fit  as  illustrated  in  Figure  4-1. 
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4.4 


Example  of  RMAFG  Program  Output 


Data  from  a set  of  computer  runs  made  as  a part  of  the  shoulder 
model  studies  (see  Section  3.2)  were  used  to  illustrate  the  output  provided 
by  the  RMAFG  computer  program.  In  this  set  of  simulations,  the  shoulders 
were  represented  as  EULER  joints  with  the  precession  and  nutation  axes 
initially  locked  and  allowed  to  unlock  at  specified  torque  levels.  The 
coulomb  friction  and  the  unlocking  torques  were  varied  over  preselected 
ranges  and  the  unlocking  torque  for  each  axis  was  arbitrarily  set  as  110 
percent  of  the  coulomb  friction.  The  response  measures  selected  for 
determining  approximating  functions  with  the  RMAFG  using  the  coulomb  friction 
of  the  precession  and  nutation  axes  as  the  independent  variables  were:  the 

maximum  pitch  angle  of  the  right  lower  arm,  the  time  at  which  this  maximum 
pitch  occurred,  and  the  pitch  and  the  yaw  angles  of  the  right  lower  arm  at 
the  end  of  each  300  millisecond  run.  The  data  from  the  shoulder  simulations 
that  were  input  to  the  RMAFG  program  are  given  in  Table  4-3. 

The  output  from  the  RMAFG  routine  for  each  of  the  four  response 
measures  evaluated  is  presented  in  Figure  4-2 (a)  through  (d) . Note  that, 
for  each  case,  a quadratic  fit  of  the  data  using  MODE  = 2 (see  Section  4.2) 
was  specified.  Also,  the  values  of  the  functions  for  two  intermediate  values 
of  the  independent  variables  were  calculated. 

A sample  plot  of  the  pitch  angle  of  the  lower  arm  at  300  msec  as  a 
function  of  the  coulomb  friction  on  the  precession  and  nutation  axes  that 
results  from  the  quadratic  fit  obtained  with  the  RMAFG  is  depicted  in 
Figure  4-3. 
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RESULTS  OF  SHOULDER  SIMULATIONS  INPUT  TO  THE  RMAFG 
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SHLD25  RUNS  TO  DEMONSTRATE  RMAFG,  PENDULUM  TEST  14  SEGMENT  SHOULDER  ARMS, EULER 
ND.VAR.  COULOMB  TORQUE;  DEP.  MAX  PITCH,  TIME  MAX  PITCH,  PITCH(300MS) ,YAW(300MS) 


Mean 

-1. 

18750 

Sigma 

2.54957 

for  Function  Y 

J Degree 

Mean 

Sigma : 

for  Parameters  X( 

) 

1 2 

141, 87500 

27.21414 

2 2 

143.25000 

25.81061 

MEAN  removed 

from  Parameters 

X(  ) and 

normalized 

The  error  using  1 coefficients  is 

2.54957 

Coefficients 

Value 

Constant 

0. 

-1 

. 187500 

The  error  using 

3 coefficients  is 

.39268 

Coefficients 

Value 

0. 

-1.187500 

1. 

-2.527440 

2.- 

.307646 

Linear 


The  error  using 

6 coefficients  is 

. 15472 

Coefficients 

Value 

0. 

-2.098679 

1. 

-2.752310 

2. 

.256477 

11. 

.768566 

21. 

-.247936 

22. 

. 164377 

Quadratic 


Error  Distribution  - Histogram 


Sigma 

If  Points 

Point  IndeX 

■2.50 

0. 

■2.00 

0. 

■1.50 

1.  » 

2 

• 1.00 

0. 

-.50 

3.  *** 

3 

5 8 

0.00 

1.  * 

7 

.50 

1.  • 

6 

1.00 

0. 

1.50 

2.  *• 

1 

4 

2.00 

0. 

2.50 

0. 

(a) 

Function 

For  Maximum  Pitch  Angle 

Figure  4-2 

SAMPLE  OUTPUT  OF  THE 

RMAFG  ROUTINE  APPLIED 

TO 

SHOULDER 

MODEL  RESPONSE  MEASURES 
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Function  1 Evaluations  at  Given  Input  Points 


Point 

Y(Input) 

Y(Fit) 

Parameters  X(  ) 

1 

-3.910 

-3.713 

163.000 

122.000 

2 

1 1020 

.752 

122.000 

163.000 

3' 

-1.900 

-1.985 

140.000 

110.000 

4 

1.900 

2. 113 

1 10.000 

140.000 

5 

-4.250 

-4.356 

180.000 

130.000 

6 

-.010 

. 101 

130.000 

180.000 

7 

-4.250 

-4 . 242 

180.000 

180.000 

8 

1.900 

1.830 

1 10.000 

121.000 

Function  1 Evaluations  at  Specified  Input  Points 
Point  Y(Fit)  Parameters  X(  ) 

1 -.895  130.000  130.000 

2 -2.793  150.000  150.000 


Figure  4-2  (a)  (Cont'd.) 


SHLD25  RUNS  TO  DEMONSTRATE  RMAFG,  PENDULUM  TEST  14  SEGMENT  SHOULDER  ARMS, EULER 
ND. VAR,  COULOMB  TORQUE;  DEP.  MAX  PITCH,  TIME  MAX  PITCH,  PITCH(30GMS) , YAW( 300MS) 

Mean  66.25000  Sigma  20.11685  for  Function  Y 


J Degree  Mean 

1 2 141.87500 

2 2 143.25000 


Sigma:  for  Parameters  X(  ) 

27.21414 
25.81061 


MEAN  removed  from  Parameters  X(  ) and  normalized 

The  error  using  1 coefficients  is  20.11685 

Coefficients  Value 


0.  66.250000 


The  error  using 
Coefficients 

0. 

1. 

2. 


3 coefficients 
Value 

66.250000 

-19.816600 

3.955115 


is  3.17906 


The  error  using  6 coefficients  is  1.43056 

Coefficients  Value 


0. 

1. 

2. 

11. 

21. 

22. 


54.792162 
-21.930729 
2.442389 
5.207368 
-5.730120 
6. 753*147 


Error  Distribution  - Histogram 

Sigma  it  Points  Point  IndeX 


-2.50 

0. 

-2.00 

0. 

-1.50 

1.  * 

2 

-1.00 

1.  * 

3 

-.50 

1.  « 

5 

0.00 

2.  »» 

7 

.50 

0. 

1.00 

2.  »« 

4 

1.50 

1.  * 

1 

2.00 

0. 

2.50 

0. 

(b)  Function  For  Time  of  Occurrence 
of  Maximum  Pitch  Angle 


Figure  4-2  (Cont'd.) 
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Function  2 Evaluations  at  Given  Input  Points 


Point 

Y(Input) 

Y (Fit ) 

Parameters  X(  ) 

1 

45.000 

47. 135 

163.000 

122.000 

2 

85.000 

82.61 1 

122.000 

163.000 

3 

65.000 

63.881 

140.000 

1 10.000 

4 

85.000 

86.577 

1 10.000 

140.000 

5 

40.000 

38.936 

180.000 

130.000 

6 

85.000 

86.082 

130.000 

180.000 

7 

40.000 

40.028 

180.000 

180.000 

8 

85.000 

84.750 

1 10.000 

121.000 

Function  2 Evaluations  at  Specified  Input  Points 
Point  Y(Fit)  Parameters  X(  ) 

1 64.596  130.000  130.000 

2 49.362  150.000  150.000 


Figure  4-2 (b)  (Cont’d.) 


SHLD25  RUNS  TO  DEMONSTRATE  RMAFG,  PENDULUM  TEST  14  SEGMENT  SHOULDER  ARMS, EULER 
ND. VAR.  COULOMB  TORQUE;  DEP.  MAX  PITCH,  TIME  MAX  PITCH,  PITCH ( 300MS ) ,YAW(3GOMS) 


Mean 


-67.85000  Sigma 


5.97668  for  Function  Y 


J Degree  Mean 

1 2 141.87500 

2 2 143.25000 


Sigma;  for  Parameters  X(  ) 
27.21414 
25.81061 


MEAN  removed  from  Parameters  X(  ) and  normalized 


The  error  using  1 coefficients  is  5.97668 

Coefficients  Value 


0.  -67.850000 


The  error  using  3 coefficients  is  2.41994 

Coefficients  Value 


0.  -67.850000 

1.  -5.470639 

2.  .071163 


The  error  using  6 coefficients  is  .56663 

Coefficients  Value 


0. 

1. 

2. 

11. 

21. 

22. 


-66.643753 
-4.942339 
-. 132234 
-2.443365 
-.556155 
1.285935 


Error  Distribution  - Histogram 
Sigma  # Points 


Point  IndeX 


-2.50 

-2.00 

-1.50 

-1.00 

-.50 

0.00 

.50 

1.00 

1.50 

2.00 

2.50 


0. 

0. 

1.  * 
1.  * 
1.  * 
2.  *» 
1.  * 
1.  » 
0. 

1.  * 

0. 


3 
2 

5 

4 7 

6 
8 

1 

(c)  Function  For  Pitch  Angle 


@ 300  msec. 


Figure  4-2  (Cont'd.) 
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Function  3 Evaluations  at  Given  Input  Points 


Point 

Y(Input) 

Y(Fit) 

Parameters  X(  ) 

1 

-71.660 

-70.617 

163.000 

122.000 

2 

-62.820 

-63.375 

122.000 

163.000 

3 

-63.280 

-64.060 

140.000 

110.000 

4 

-64. 150 

-64.252 

1 10.000 

140.000 

5 

-77. 140 

-77.556 

180.000 

130.000 

6 

-62.570 

-62. 188 

130.000 

180.000 

7 

-77.000 

-77.054 

180.000 

180.000 

8 

-64. 180 

-63.699 

110.000 

121.000 

Function  3 Evaluations  at  Specified  Input  Points 
Point  Y(Fit)  Parameters  X(  ) 

1 -64.670  130.000  130.000 

2 -68.327  150.000  150.000 


Figure  4-2 (c)  (Cont'd.) 
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SHLD25  RUNS  TO  DEMONSTRATE  RMAFG,  PENDULUM  TEST  14  SEGMENT  SHOULDER  ARMS, EULER 
ND. VAR.  COULOMB  TORQUE;  DEP.  MAX  PITCH,  TIME  MAX  PITCH,  PITCH(300MS) , YAW( 300MS) 


Mean 

10. 10500 

Sigma 

2.80935 

J 

Degree 

Mean 

Sigma: 

for  Parameters  X(  ) 

1 

2 

141.87500 

27.21414 

2 

2 

143.25000 

25.81061 

MEAN 

removed 

from  Parameters 

X(  ) and 

normalized 

The 

error  using  1 coefficients  is 

2.80935 

Coefficients 

Value 

0. 

10. 105000 

The  error  using 

3 coefficients  is 

.98164 

Coefficients 

Value 

0. 

10.  105000 

1. 

2.603300 

2. 

.223047 

The  error 

using 

6 coefficients  is 

. 10535 

Coefficients 

Value 

0. 

7.943752 

1. 

2.052303 

2. 

.070976 

11. 

1.855552 

21. 

-.310125 

22. 

.332918 

Error  Distribution  - Histogram 

Sigma 

It 

Points 

Point  IndeX 

-2.50 

0. 

-2.00 

0. 

-1.50 

1. 

» 

2 

-1.00 

0. 

-.50 

3. 

*«« 

3 5 8 

0.00 

1. 

» 

7 

.50 

1. 

» 

6 

1.00 

1. 

* 

1 

1.50 

1. 

« 

4 

2.00 

0. 

2.50 

0. 

(d)  Function 

For  Yaw  Angle 

@ 300  msec. 

Figure  4-2  (Cont'd.) 
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Function  4 Evaluations  at  Given  Input  Points 


Point 

Y(Input) 

Y(Fit) 

Parameters  X(  ) 

1 

10.910 

11.020 

163.000 

122.000 

2 

8.040 

7.857 

122.000 

163.000 

3 

8.280 

8.245 

140.000 

1 10.000 

4 

7.870 

8.036 

110.000 

140.000 

5 

14.800 

14.735 

180.000 

130.000 

6 

8.300 

8.370 

130.000 

180.000 

7 

14.610 

14.618 

180.000 

180.000 

8 

8.030 

7.959 

110.000 

121.000 

Function  4 Evaluations  at  Specified  Input  Points 
Point  Y(Fit)  Parameters  X(  ) 

1 7.383  130.000  130.000 

2 8.739  150.000  150.000 


Figure  4-2 (d)  (Cont'd.) 
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Figure  4-3  SAMPLE  PLOT  OF  RMAFG  QUADRATIC  FUNCTION 
FOR  PITCH  ANGLE  @ 300  MSEC. 


SO 


5.0 


REDUCTION  OF  THE  NUMBER  OF  INERTIA  RELATED  PARAMETERS 


A fact  that  is  not  widely  appreciated  is  that  there  is  an  infinite 
number  of  different  chains  of  connected  rigid  bodies,  each  with  different 
mass  distributions,  that  all  have  identical  dynamic  responses  and  may  be 
considered  as  "dynamically  equivalent."  In  particular,  it  is  possible  to 
add  or  subtract  an  arbitrary  amount  from  each  mass  in  a connected  chain 
of  rigid  bodies  and  obtain  a dynamically  equivalent  system.  The  only  con- 
straint is  that  the  sum  of  the  mass  perturbations  be  zero,  i.e.,  the  total 
mass  must  not  be  changed.  It  is,  in  fact,  possible  to  make  one  or  more  of 
the  masses  in  the  chain  zero,  or  negative,  and  still  define  an  equivalent 
system. 


The  purpose  of  this  task  was  to  develop  a computer  program  to 
compute  dynamically  equivalent  systems.  The  basic  theorem  and  the  questions 
are  given  on  pages  113-117  of  Volume  1 - Engineering  Manual  of  Reference  1. 

Theorem: 


Given  a system  of  rigid  bodies  connected  by  joints  in  a tree 
structure,  make  the  following  transformations: 


M,  = where  Z = 0 mass  equation, 

k=l 


pk 

" Fk  + Z.  dk. . 

x (v^  x 

],  inertia  equation. 

* 

xk 

1 

II 

'j  J , 

transformation 

of  c.g. , 

Mkck 

= Z d.  r, 

. k . k . 

t- 

evaluation  of 

cv 

* 

rk. 

• J J J 

= rk.  + ck> 

transformat  ion 

of  joint  locations. 

J J 


(The  notation  [s  x (r  x ] denotes  the  matrix  vs  - s vl3 

T 

where  s is  the  transpose  of  s and  I is  the  identity 
matrix. ) 
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where : 


"k 


dk 

dk 


P 


k 


.+ 

is  the  mass  of  segment  k and  My  is  the 
perturbed  mass, 

is  the  mass  perturbation  of  segment  k3 

is  the  sum  of  the  dy  of  the  segments  which  may 
be  reached  through  joint  j from  segment  k3 

* 

is  the  inertia  matrix  of  segment  k.  and  Py 
is  the  perturbed  inertia. 


X y is  the  location  of  the  c.g.  of  segment  k and 

* 

Xy  is  the  perturbed  location, 

C y is  the  perturbation  of  the  c.g.  and  joint  locations 

of  segment  k3 


v y is  the  location  of  the  j'th  joi.it  associated  with 

J segment  k and 

* 

v ^ is  the  perturbed  joint  location.  Joint  locations 

are  relative  to  the  c.g.  of  the  segment. 


The  theorem  states  that  this  perturbed  system  has  the  same  motion  as 
the  unperturbed  system.  This  implies  that,  even  if  one  had  complete  knowledge 
of  the  motion,  i.e.,  positions,  velocities  and  orientations  of  the  segments,  one 
would  not  be  able  to  compute  the  masses  of  the  individual  segments  (the  My)  or 
determine  the  location  of  the  center  of  mass  of  the  segments  (the  J,) . Note 

K 

that  the  geometry  is  the  same.  The  mass,  inertia  matrix  and  location  of  the 


center  of  mass  are  the  only  quantities  that  are  changed.  (Of  course,  the 
location  of  a joint  relative  to  the  c will  change  because  the  j*  will 
change,  but  the  distance  between  joints  attached  to  the  same  segment  is  not 
changed. ) 

5. 1 Description  of  the  Delta  Algorithm  Fortran  Program 

The  listing  of  the  Fortran  program  called  the  "Delta  Algorithm" 
which  will  compute  an  equivalent  system  is  given  in  Appendix  C.  Figure  5-1 
is  a flow  diagram  of  the  main  program. 

The  routine  to  compute  an  equivalent  system  has  been  written  so 
that  it  accepts  the  segment  cards  (B2)  and  the  joint  cards  (B3)  with  the 
same  information  and  format  as  the  main  CVS  program  (see  Reference  1,  Volume  3). 
A lead  card  with  the  number  of  segments  (NSEG) , the  number  of  joints  (NJNT) 
and  the  value  of  G (to  convert  weight  to  mass  units)  is  required.  A final 

, I 

card(s)  containing  the  weight  (mass)  perturbations  is  also  read. 

All  of  the  computations  are  made  in  subroutine  DELTA.  The  key  to 

the  procedure  is  the  computation  of  the  from  i;he  information  in  JNT(J). 

J 

The  logic  of  the  computation  is  in  the  DO  35  loop  of  the  subroutine.  Although 
it  has  only  about  25  Fortran  statements,  this  logic  is  quite  involved  and 
hence  is  difficult  to  follow.  It  is  based  on  the  following  properties  of  the 
JNT  vector. 

0 < JNT(J)  <_  J 

(This  assumes  all  segments  are  connected  and  any  flexible  elements  are 
considered  as  regular  segments.) 
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START 

I 

I 


INPUT-  Title,  segment  weights,  inertia, 
joint  array  and  locations. 

Program  accepts  cards  in  standard  CVS  format 
CVS  cards  A1  through  B3 


I 

I 


INPUT-  Print  option  and  mass  perturbations 
are  inputted  from  the  terminal. 


J PRINT-  Inputs  are  printed 


I 

I 

I 

I 


J CALL  DELTA  ALGORITHM 


I 

I 


{ CALL  EIGEN  - process  inertia  matrices 


I 

I 


J PRINT  perturbed  masses,  inertia  matrices, 

! principal  components  of  inertia  and  orientation 
! of  the  principal  axes,  new  joint  locations. 


a 

i 

STOP 


Figure  5-1  FLOW  DIAGRAM  OF  DELTA  ALGORITHM  COMPUTER  PROGRAM 
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If  JNT {J ) = N,  the  segment  J+l  is  in  the  branch  emanating  from 
segment  N via  joint  J.  All  segment  numbers  which  are  not  in  the  vector 
JNT(J)  are  the  end  of  a branch  (except  for  segment  1).  For  example,  consider 
the  JNT  vector  of  the  standard  15  segment  crash  victim  model: 

J 1 2 3 4 5 6 7 8 9 10  11  12  13  14 

JNT  123416719  10  3 12  3 14 

Segments  5,  8,  11,  13,  15  are  not  in  the  vector.  These  are  the  head, 
feet  and  the  lower  arms. 

In  this  example,  for  segment  1 


d2  + d3  + d4  + dS  + d12  + d13  + d14  + d15 


d1  = d6+d7+d8 

u 


d9  + d10  + dll 


Note : 


d~  + d-  + d-  + d^  = 0,  since  id  = 0, 
1 21  1 5 28 


In 


the  routine  the  sums  required  to  compute  and  P^  are  accumulated  as 


a^3  and  the  d ^ are  determined  in  the  DO  35  loop.  The  formula  for  P^ 


is  written  in  a modified  form  as  follows: 


* 

Pj^  — P -h  l k.  P ^ ^k  ^~^k  ^ ^k  ^ ^ 

>7  .7  .7  J 3 J J 

* 

= P.^  + l d-^  [r^,  3?  ^ ~^k  ^k  ^~Pk  ^ Pk  ^ ^ 

J 3 J 3 

T T 

(The  notations  [s  x (r  r]  denotes  the  matrix  vs-  s vl3 
T 

where  s is  the  transpose  of  s and  I is  the  identity 
matrix. ) 

* 

The  Vj^  are  computed  in  subroutine  DELTA.  When  subroutine  DELTA 

3 

has  finished,  the  main  program  calls  subroutine  EIGEN  to  determine  the 

* 

principal  moments  of  inertia  from  P ^ and  the  relative  direction  cosine 

matrix.  The  main  routine  then  computes  the  yaw,  pitch  and  roll  angles  of  the 
* 

system  for  P with  respect  to  the  original  system  for  P^. 

* A 

Finally,  the  Mk  * Pk.  J Ck  J principal  moments  of  inertia,  yaw 

pitch  and  roll,  and  the  r ^ * are  printed.  This  is  sufficient  information 

3 

to  form  an  input  deck  for  the  latest  version  of  the  CVS  program  which  allows 
for  a rotated  inertia  matrix. 


5.2  Sample  Application  if  the  Delta  Algorithm  Program 


The  Delta  algorithm  computer  program  was  executed  to  illustrate  how 
it  can  be  used  to  compute  dynamically  equivalent  systems.  For  this  example, 
an  old  data  deck  for  a 15  segment-14  joint  model  of  a crash  victim  was  used. 

The  printout  from  the  routine,  presented  in  Figure  5-2,  shows  the  values  of 
the  original  physical  system  and  the  weight  perturbations  of  the  various  segments 
that  were  input  to  the  program,  and  the  computed  characteristics  of  the  segments 
and  joints  of  the  dynamically  equivalent  system. 
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Figure  5-2  SAMPLE  OUTPUT  FROM  THE  DELTA  ALGORITHM  COMPUTER  PROGRAM 


COM PUT AT I OHS  OF  DELTA  ALGORITHM 
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Figure  5-2  (Cont'd.) 
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iCGMEMT  9 WEIGHT  17.190 

INERTIA  MATRIX  OFFSET  EIGENVALUES  DIRECTION  COSINE  Y-P-R 
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0.00000  0.00000  .19400  .90PR7  .15400  0.000000  0.000000  1.000000  0.00000 
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6.0 


ANSWERS  TO  SPECIFIC  QUESTIONS 


Task  8 of  the  research  effort  was  to  answer  a series  of  questions 
designed  to  address  the  problem  of  determining  the  required  level  of 
model  detail  (level  of  sophistication)  that  would  be  adequate  to  simulate 
typical  crash  situations  and  the  problems  related  to  the  development  of 
techniques  for  performing  simulation  studies.  The  principal  questions  were 
listed  earlier  in  Table  1-1  and  the  levels  of  sophistication  in 
Table  1-2. 


The  following  subsections  describe  the  analytical  studies  that  were 
performed  to  answer  several  of  the  questions  of  Table  1-1.  It  was 
originally  intended  to  supplement  these  studies  with  sample  runs  of  the  CVS 
model  but  these  runs  were  not  completed. 

6. 1 Table  1-1,  Question  1 

Question  1 in  Table  1-1  is:  Is  similitude  between  rigid  body 

characteristics  of  selected  segments  adequate  or  must  similitude  of  one  or 
more  vibration  modes  also  be  preserved? 

Our  answer  to  this  question  is  a qualified  yes;  rigid  body 
characteristics  are  adequate  for  all  segments  except  the  lumbar  spine  and 
neck.  The  spine  and  neck  are  discussed  elsewhere  in  this  report.  The  context 
in  which  we  are  answering  this  question  is  one  where  the  concern  is  with  the 
gross  motion  of  the  dummy  or  the  human  and  not  with  more  sophisticated  details 
such  as  the  effects  on  internal  organs  like  the  brain  or  those  contained  in  the 
torso.  Thus,  we  interpret  this  question  as  one  concerned  with  the  long  bone 
segments  which  are  the  arms  and  legs  and  answer  the  question  accordingly. 

To  substantiate  our  answer  we  rely  on  the  following  development 
which  is  based  on  the  simple  beam  theory  given  in  Reference  3. 
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Consider  a simply  supported  beam  (one  in  which  the  ends  are  held  by 
pin  or  ball  joints)  as  depicted  in  Figure  6-1.  The  frequencies  of  free 
oscillation  of  such  a beam  are  given  by  (Reference  3,  page  83): 


2 

n it 


where  = frequency,  Hz 


= mode  number  (1,2,----) 


length  of  beam,  in. 


E 

I 

P 

A 

G 


= Young's  modulus  (30x10  lb/in  for  steel) 

4 

= second  moment  of  area,  in 

= density  (0.28  lb/in2 3  for  steel) 

-2 

= cross  sectional  area,  m 

2 

= 386  in/sec 


PHt) 


Y 


Figure  6-1  SIMPLY  SUPPORTED  BEAM 
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For  a circular  rod  I = 0 . 5 tt  r and  A = tt  r where  v is  the 
radius.  Hence, 

f = 2.25  x 105  n r/l2. 

J n 

2 

If  l =12  in.  and  r = 1/4  in.,  f = 390  n . Thus,  the  fundamental 

J n 

frequency  of  oscillation  (n  = 1)  is  390  Hertz,  or  a 2.6  millisecond 
period.  The  first  harmonic  (n  = 2)  would  be  1560  Hertz  or  a 0.64  millisecond 
period.  The  mode  shape  (displacement  of  rod)  is  sinusoidal; 

<f>  (x)  = sin  ( nirx/H ) , 0 <_  x <_  l . 

If  the  beam  is  subjected  to  an  axial  compressive  load  of  Q lbs. 
(tensile  load  -Q) , the  frequency  is  reduced  (increased)  by  the  factor 

-y/  1-Q/Q  where  Q is  the  Euler  critical  load  (Reference  3,  page  136).  The 

° ° 2 2 

Euler  critical  load  is  given  by  = tt  EI/l  . For  a steel  circular  rod  12 

inches  long  and  0.25  inch  radius;  Q = 12600  lb. 

a 

Response  to  Disturbing  Force 

If  the  beam  is  subjected  to  a normal  disturbing  force  located  at 
one  point  (see  Figure  6-1)  the  displacement  y is  given  by: 

u = Z <p  (x)  q (t) 

^ n ^ n 

n 

where : 

••  2 

q + a)  q = P <p  (a)  f(t)/m 
^ n n rn 


(x) 

n 


sin 


fmrx/ijj  mode  shape 


P f(t)y  disturbing  force 


m = p Ai , mass  of  beam. 
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If  the  disturbing  force  is  sinusoidal,  f(t)  = sin  ut , the  solution 

2 2 

is  n = a sin  w t + b cos  a > t P <p  (a)  sin  ut/(m(u  - w ) ) where  a 
in  n n n n n n n 

and  b are  constants  depending  on  the  initial  conditions. 

The  steady  state  response  will  be 

y = y(x)  sin  ut  where 

y(x)  = (P/M)l  <p  (a)  <p  fxj/fui2  - u2) . 
n 

2 2 

When  the  frequency  is  small  (w  <<u^) j y(x)  will  approach  the 
deflection  due  to  a static  force  P at  a as  given  by: 

y(x)  = (Pi2 /El)  h(x)  (6.1) 

where : 

2 2 4 

h(x)  = (i-a)  x (2ia  - a -x)/(6i)f0<_x<_a 

= (i-x)  a (i2  - o'  - (i  - x)2 ) / ( 6i4 ) , a <_  x <_  i 
The  maximum  static  deflection  occurs  at  = a and  is 
y(a)  = (Pi2 /El)  (a.2 (i-a)2 /Zi4 ) 

= (Pi2  )/4  8 El  if  a = i/2 

For  example,  for  a circular  steel  rod  12  inches  long  and  0.25  inch 
radius,  if  the  force  is  applied  at  the  mid  point  the  static  deflection  is 
0.2  inches  per  1000  lbs  of  force. 
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Force  and  Moment 


The  moment  M acting  at  the  point  x is: 
M =-EI 


and  the  shear  force  S is : 


S 


3 M 
3x 


From  the  equation  6.1,  the  moment  M is  zero  at  the  end  points  and  equal  to 
P(l-a)  at  x = a.  The  shear  force  S is  P(l-a)/l  from  x = 0 to  x = a; 

and  -Pa/l  from  x = a to  x = l.  The  discontinuity  at  x = a is  due  to  the 

applied  load  P at  that  point.  For  the  case  of  free  oscillations  where  the 

mode  shape  is  sinusoidal  we  have  for  a unit  peak  displacement: 

„ 2 2 , 

pt  3 • /ttTrx,  ,m\  > . ,m\x/ . , , 

M = -PI  — - s^7'2  (— — ) = El  (—)  s%n  ( l)  and 

A U X/  X/ 

9x 

3 M „T  ,m\  ,3  ,nvx/ . , 3n  t r ,mrx/n. 

S = — = El  (—)  aos  ( l)  = n Q — jos  ( i) 

3 x l cl 

where  Q is  the  Euler  critical  load.  At  x = 0 for  a circular  steel  rod 

c 

12  inches  long  and  0.25  inch  radius,  S = 3300  lb  per  inch  peak  displacement 
of  the  first  mode  (n  = 1) . 

Discussion 


The  preceding  development  was  concerned  with  effects  of  transverse 
loading.  Axial  loads  will  be  propagated  with  the  speed  of  sound.  In  steel 
the  speed  of  sound  is  2 x 10^  inches  per  second.  This  is  60  microseconds  per 
foot.  Since  the  typical  minimum  time  interval  of  interest  in  the  gross  motion 
simulation  is  the  order  of  1 millisecond  and  typical  dimensions  are  of  the 


order  of  1 foot,  propagation  time  delays  can  reasonably  be  ignored.  The 
compression  or  expansion  due  to  axial  loads  will  be  negligible  for  loads 

experienced  in  typical  situations.  The  reduction  in  length  due  to  bending 

2 

from  a transverse  load  will  be  of  the  order  of  d /£  where  d is  the  peak 
displacement  and  £ is  the  length. 

The  transverse  loading  considered  assumed  a simply  supported  beam 
which  implies  that  the  ends  of  the  beam  are  held  by  a pin  or  a ball  to 
adjoining  segments  which  are  held  fixed.  In  the  typical  situation  this  is 
not  true  but  they  may  be  partially  'fixed'  by  contact  loading  and/or  inertial 
loading.  Any  of  these  situations  will  reduce  the  peak  displacement.  These 
may  be  estimated  from  the  shear  forces  which  exist  at  the  ends  of  the  beam. 

Conclusion 

For  the  long  bone  segments  in  the  Part  572  dummy  the  deformation 
due  to  axial  loading  is  negligible  and  the  deformation  due  to  transverse 
loading  will  probably  not  exceed  a small  fraction  of  an  inch  per  1000  pounds 
of  load.  It  is  believed  that  errors  caused  by  ignoring  these  effects  will 
be  much  smaller  than  the  errors  caused  by  other  simplifying  assumptions  in 
the  model. 

6.2  Table  1-1,  Questions  5 and  7 

Questions  5 and  7 of  Table  1-1  concerning  required  level  of  model 
detail  are: 

5:  Is  it  sufficient  to  preserve  similitude  of  impulses  and 

coefficients  of  restitution  during  'hard'  impact  processes  or  must  simulitude 
of  forces  also  be  preserved? 

7:  Is  a sliding/rolling  characterization  of  a particular  contact 

adequate  or  must  similitude  of  the  deflection  characteristics  be  preserved? 
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To  answer  these  questions  consider  the  idealized  problem  of  a hard 
sphere  contacting  a deformable  planar  surface  where  the  force-deflection 
characteristic  of  the  contact  is  linear  and  the  coefficient  of  friction  is  a 
constant.  The  geometry  of  the  contact  is  illustrated  in  Figure  6-2  below. 

Z 


(A7/1X/A7UM  P£NET/?A  T/OAJ ) 
Figure  6-2  SPHERE-PLANE  CONTACT 


Let  the  motion  be  jn  the  X-Z  plane.  The  equations  of  motion  are: 


MX 

— 

-pk  (r-z) 

(6.2) 

MZ 

= 

k (r-z) 

Mp^Q 



pk.r  'r-z) 
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where 


where 


M 

- 

mass 

lb-sec^/in 

P 

= 

radius  of  gyration 

in 

X 

= 

X position  of  center  of  sphere 

in 

Z 

= 

Z position  of  center  of  sphere 

in 

x/z 

— 

linear  velocities 

in/sec 

2 

Z,Z 

= 

linear  accelerations 

in/sec 

0 

= 

angular  position  of  sphere 

radians 

0 

= 

angular  velocity 

rad/ sec 

• • 

0 

= 

angular  acceleration 

rad/sec^ 

V 

= 

radius  of  sphere 

in 

k 

= 

spring  constant 

lb/in 

y 

- 

friction  coefficient 

The  solutions  of  equations  6.2  are: 

t = 0 


(6.3) 


ti  = \k/M 


X = Xq  + vZq  (1  - oos  (tit)) 


initially  at 


Z = Z^  oos  tit  s 
0 = 0^-  yrZ^  (1  - oos  (tit))/p^  j 
X = Xgt  + \iZq  (t  - sin  (tit) /ti)  j 
Z = r + Zq  sin  tit/ti  3 
0 = Q t - yrZ^  (t  - sin(tit) /ti) /p2  j 


0 


r 

0 
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At  time  equal  to  zero,  the  time  of  initial  contact,  assume  that 
Z^  is  negative  and  Xp  is  positive.  The  motion  in  the  Z direction  is  not 
affected  by  the  friction.  At  time  t = t\/2SI,  Z = 0 and  the  value  of  Z is 
v + Z^/ft;  this  is  the  time  of  maximum  penetration  which  is  Z^/ft  giving  a 
maximum  spring  force  of  k(r-z)  = At  time  t = -rr/fi,  the  sphere  has 


rebounded  and  leaves  the  surface  with  a Z velocity  of  Z = -ZQ.  This  is 


analogous  to  an  impulsive  impact  with  a coefficient  of  restitution  of  unity. 


The  behavior  in  the  X direction  is  dependent  on  the  coefficient  of 
friction.  The  tangential  ( X direction)  velocity,  V,  at  the  point  of  contact 
is  given  by 


2,  2 


V = XQ  - rdg  + \i  (1  -hr/ p J Z^  (1  - cos  (tit)) 


(6.4) 


The  sign  of  the  friction  coefficient  \i  is  selected  to  agree  with  the  sign  of 
V,  i.e.,  the  frictional  force  opposes  the  tangential  motion.  Initially 


VQ  = XQ  - rQg.  We  must  distinguish  between  several  cases  which  are: 


Case  1: 


The  tangential  velocity,  V,  is  initially  positive  and  remains 
positive  throughout  the  interval  0 <_  t <_  tt /ft  3 T. 


Case  2: 


' 

I 


The  tangential  velocity,  Va  is  initially  positive  and  becomes 
zero  at  time  t = vT  where  0 < u < 1 and  u T satisfies  the  equation 


V = XQ  - rdQ  + y (1  + r2/p2)  ZQ  (1  - cos  (mT))  = 0. 


Case  3: 


The  tangential  velocity,  V3  is  initially  zero.  In  this  case  (and 

also  in  Cases  2 and  5 when  the  tangential  velocity  becomes  zero) 

• • •» 

the  tangential  velocity  remains  zero  since  X = 0 = 0 and  the 
sphere  rolls  on  the  plane. 
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Case  4: 


Case  5 


leaves 


where 


and 


The  tangential  velocity,  V3  is  initially  negative  and  remains 
negative  throughout  the  interval  0 <_  t <_  T . In  this  case  and  in 
Case  5,  -y  is  substituted  for  y in  equations  6.2,  6.3  and  6.4. 

The  tangential  velocity,  V3  is  initially  negative  and  becomes 
zero  at  time  t = u T3  where  0 <_  u <_  1 and  u T satisfies  the 
equation 

V = XQ  - rhQ  - \x  (1  + r2/p2)  ZQ  (1  - cos  (vT) ) = 0. 

The  solution  to  all  five  cases  at  t = T3  the  time  when  the  sphere 
the  surface,  can  be  written  as: 

X = Xq  + \iZq  (1  - cos  (wn)) 


0=0^-  yrZ^  (1  - cos  (vu))/p2  (6.5) 

X = (Xq  + yZ q)  T - V-ZqT  (sin  (vttJ/tt  + (1- o)  cos  fu it)) 

Z = r 


• % 9 • 9 

0 = (Qq  -nrZg/p  )T  f yrZ^T  (sin  (vttJ/tt  + (1-v)  cos  (wn))/p 


y is  replaced  by  -y  if  vo  • <xo  - rV  < 


V=0  if  vfl  - I - P90  - o. 


v is  the  solution  of  the  equation: 


V = XQ  - rQQ  + p (1  + r2/p2)  ZQ  (1  - cos  (vtt))  = 0 


(6.6) 
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providing  the  solution  is  in  the  range  0 <_  u £ 1, 


u = 0 if  Xq  - vQq  = 0, 


u = 1 if  ^ + rQ0  ^ 0 and  equation  6.6  does  not  have  a solution 
in  the  range  0 <_  u <_  1. 

Impulsive  Solutions 

Consider  the  case  of  the  sphere  contacting  the  plane  where  we  apply 
an  impulse  at  the  time  of  first  contact.  The  equations  are: 

MLX  = -pi 

MAZ  = J 

2 • 

Mp  A6  = prl 

• • • 

where  AJ,  AZ,  and  A9  are  the  instantaneous  changes  in  velocity  due  to  the 
impulse  of  strength  I (lb-sec). 


impulse) 


Letting  I be  equal  to  -2b IZq  we  have  a.:  t = 0 + (just  after  the 


z+  - -zo 


Q+=  Q Q - 2prZQ/p‘ 
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where  the  sign  of  the  friction  coefficient  is  the  same  as  the 
relative  velocity  V ^ = X ^ - rQ ^ at  the  time  of  contact.  If 
taken  as  equal  to  zero.  At  a time  T we  have 


sign  of  the 

VQ  = 0,  y 


is 


* - *0  + 2vi0 

z = -ig 

0 = - 2urZg/p2  (6.7) 

X = (XQ  + 2\iZ  )T 
Z = v -ZqT 

e = (Q0  - 2]irZ0/p2)T 

Comparison  of  Solutions 

The  results  of  the  'soft'  contact  are  given  by  Equations  6.5  and  of 
the  'hard'  contact  byEquations  6.7  at  time  T when  the  'soft*  contact  has  left 

the  surface.  For  the  case  of  pure  roll  (Case  3),  V.  = 03  we  have: 

!/ 

'soft*  ’hard'  (impulse) 


X 

xo 

*0 

• 

• 

z 

-zo 

• 

* 

• 

0 

% 

X 

V 

V 

z 

r 

r-Z0T 

• 

e 

CD 

‘-3 

V 
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We  note  that  the  only  difference  is  in 
'soft'  contact  becomes  harder  (k  increases),  T 
T ■+■  0 and  the  results  are  identical.  For  Cases 


the  Z position.  As  the 
gets  smaller.  In  the  limit, 

1 and  4 where  u = 1 we  have 


' soft  * 


' hard ' 


X 

z 

e 

x 

z 

0 


xo  + 2»zo 

-L 


Qq  - 2\srZQ/ p‘ 
(X0  + uz  0)T 


r 

(QQ  - \irZ0/p2)T 


XQ  -f  2uZ^ 

-*o 

eQ  - 2\irZQ/p2 
(XQ  + 2pZ0)T 
r - ZqT 

fbQ  - 2urZQ/p2)T 


We  note  that  the  velocities  are  the  same  but  the  positions  are 
different.  The  differences  in  X and  0 vanish  as  y 0 or  as  T 0 . 
The  difference  in  Z vanishes  only  when  T ■+  0. 


In  Cases  2 and  5 where  the  slide  becomes  a roll  at  some  point  in  the 
time  interval  T (0  < u < 1)  there  will  be  a difference  in  both  the  velocities 
and  positions.  For  example  for  u =1/2 

'soft'  'hard' 


X 

z 

0 

X 

z 

0 


0 - yrVp2 

ri0  + VZC)T  - VigT/v 
V 

(Qq  - \irZQ/p2 ) T 
+ \irZgT/(  up2) 


0O  - 2prZ/p2 

(XQ  + 2\iZ0)T 
r - ZqT 

(bQ  - 2prZQ/p2)T 
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Numerical  Considerations 


The  time  interval  T is  given  by  T = m/ft  and  & = k./m  . 

The  mass  m is  the  weight,  w3  in  pounds  divided  by  the  acceleration  of 

2 

gravity  g (in/sea  ).  Values  of  T and  ft  for  a range  of  stiffness  to 
weight  ratio  are  given  below: 


k/w  (1/in) 


T (milliseconds,) 


ft  (rad/sec) 


1 

160 

10 

50.6 

100 

16 

1000 

5.1 

10000 

1.6 

19.6 

62 

196 

621 

1960 


Consider,  for  example,  an  impact  with  k/w  = 1000.  In  this  case  the 
sphere  would  lose  contact  with  the  planar  surface  at  T = 5.1  milliseconds.  If, 
on  the  other  hand,  an  impulsive  contact  was  assumed  and  the  impact  occurred 
at  a velocity  of  30  mph  (Z  = 528  in/sec) , the  Z position  of  the  sphere  at 
that  time  would  be  528  (.0051)  = 2.7  inches  above  the  planar  surface. 

Integrator  Considerations 


For  the  case  of  a 'soft'  impact,  the  time  T represents  a half  cycle 
of  a sinusoid.  If  we  assume  that  at  least  10  integration  steps  are  needed  to 
achieve  a reasonable  precision  in  the  numerical  integration,  then  the  maximum 
step  size  of  the  integrator  must  be  T/ 10.  The  use  of  a 'hard'  impact  imposes 
no  restrictions  on  the  integration  step  size. 

Injury  Criteria 

A commonly  used  injury  criterion  is  the  HIC  number  which  is  defined 
by  the  formula 
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2.5 


HIC 


max 


t2itl 


(t2-t2) 


1 


a 


(6.8) 


where  \a  I is  the  magnitude  of  the  acceleration  expressed  in  g's.  For  the 
'soft'  contact  we  have 

|a|  = Z0  | ^ 1 + y2  sin  (Q.t)/g 


The  HIC  number  is  computed  for  = 0.518  and  = T - From 

Equation  6.8  we  get 


HIC  = 1.3  ( 


1 + 


2 . ,2. 

y /g) 


5 n 1.5 


For  the  'hard'  impact  the  HIC  number  is  not  defined  since  by 
definition  of  an  impulse  the  integral  in  Equation  6.8  is  finite  and  the  HIC 
varies  inversely  as  (£  - Thus  the  HIC  number  would  be  infinite. 


Another  commonly  used  criterion  is  the  Severity  Index  which  is 

defined  as 


SI 


For  the 
SI  = ( 


soft'  impact  this  is  given  by 

/IT 

, , . ,2.5 

a\i  ( s%n  ]i) 

0 


where 


0 


d\i 


(sin  \x) 


2.5 


1.44 
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Thus,  for  a half  sine  wave  acceleration  the  HIC  and  the  SI  are  essentially 
equal.  For  a 'hard'  impact  the  SI  is  infinite. 

For  a specified  HIC  number  we  can  compute  fi  as  a function  of 


Z^  and  y.  The  results  for  a HIC  = 


ZQ  (in/sec) 

2 qT  (in) 

50 

0.062 

100 

0.393 

150 

1.161 

200 

2.500 

250 

4.534 

300 

7.370 

1000  and  y = 0 

are  shown  below 

n (rad/sec) 

k/w  (1/m) 

2532 

16612 

799 

1653 

406 

427 

251 

163 

173 

78 

128 

42 

Conclusions  and  Answers  to  Questions  5 and  7 of  Table  1-1 

The  only  definite  answers  we  can  make  to  the  questions  is  that  if 
the  contact  involves  a segment  where  the  HIC  number  or  the  severity  index  is 
required  the  impulse  option  or  the  slide/roll  option  must  not  be  used. 

(Proper  use  of  the  slide/roll  option  requires  the  use  of  an  impulse  to  reduce 
the  normal  velocity  to  zero  at  the  time  of  first  contact,  hence  the  slide/ 
roll  option  is  subject  to  the  same  restrictions  as  the  impulse  option.)  In 
contacts  where  the  HIC  number  or  the  severity  index  are  not  important  for 
the  segment  involved  in  the  contact,  the  user  must  make  the  decision  based  on 
evaluation  of  the  differences,  primarily  in  positions,  of  the  variables 
involved  and  the  computing  time.  It  should  be  remembered  that  the  principal 
advantage  of  the  impulse  options  in  the  CVS  program  is  that  they  place  no 
restriction  on  the  step  size  of  the  integrator  whereas  the  use  of  a 'soft* 
impact  restricts  the  step  siz.e. 
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6.3 


Table  1-1,  Question  9 


Question  9 of  those  listed  in  Table  1-1  is: 

Must  the  similitude  of  surface  compliances  be  preserved  in  the 
region  of  belt  or  air  bag  contacts? 

The  answer  to  this  question  is  yes.  To  support  this  answer  we  refer 
to  the  analytic  studies  presented  in  subsections  6.3.1  and  6.3.2  which  follow. 

A simple  model  of  a belt  passing  over  a deformable  spherical  surface 
is  analyzed  in  Section  6.3.1.  In  this  analysis  it  is  shown  that  the  effective 
stress  is  significantly  reduced  from  what  would  be  computed  if  the  compliance 
of  the  surface  was  ignored.  When  using  the  belt  algorithm  in  the  CVS  program 
the  user  must  define  a stress-strain  function  which  accounts  for  the  effect  of 
the  compliance. 

The  analysis  also  solves  the  equations  for  a belt  on  a nondeformable 
spherical  surface  with  friction,  and  shows  that  the  belt  will  not  necessarily 
lie  in  a plane  if  there  is  friction  perpendicular  to  the  belt  line.  It  also 
shows  that  the  tension  decreases  exponentially  along  the  belt  line.  The  belt 
algorithm  in  the  CVS  program  should  not  be  used  if  the  effects  of  friction  are 
important  (this  algorithm  assumes  either  zero  or  infinite  friction).  However, 
the  harness  algorithm  allows  for  surface  compliance  and  finite  friction. 

For  air  bag  restraints,  the  nominal  pressure  in  the  uncontacted  bag 
is  usually  a few  inches  of  water  and  contacts  by  body  segments  reduce  the 
volume  of  the  bag  so  that  pressures  of  a few  pounds  per  square  inch  are  pro- 
duced. The  compliance  of  the  bag  must  be  considered  to  predict  these 
pressures.  Also,  the  effecti.ve  area  is  larger  than  the  actual  area  of  contact 
as  illustrated  in  Figure  6-3. 


Figure  6-3  EFFECTIVE  AREA  OF  SEGMENT  CONTACT  WITH  AIR  BAG 


The  air  bag  algorithm  in  the  program  estimates  the  reduction  in 
volume  and  the  effective  area  assuming  that  the  contacting  segment  is  rigid 
and  that  the  bag  material  does  not  stretch.  Actual  experiments  with  an  air 
bag  (Reference  1 - Volume  2)  indicated  that  although  the  fibers  in  the  cloth 
of  the  bag  do  not  stretch  the  material  effectively  stretches  at  an  angle  to 
the  weave  so  as  to  produce  undulations  of  the  bag  surface  in  the  vicinity  of 
the  contact  by  an  object.  The  analyses  described  in  Section  6.3.2  show  this 
effect.  One  of  the  studies  shows  how  a fabric  with  a rectangular  weave  will 
lie  on  a sphere.  The  other  study  shows  the  effect  of  bag  stretch. 

6.3.1  Belt  Analyses 

(a)  Belt  With  Friction 


sketch. 


Consider  a belt  lying  on  a surface  as  shown  in  the  following 
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let : 


Assume  the  belt  may  be  treated  as  a curve  in  space.  At  the  point 


T 

N 

n = 

<7 

f 

u = 

ds  = 

4> 


unit  vector  tangent  to  belt  (in  direction 
of  decreasing  tension) 

unit  vector  principal  normal  to  belt 

unit  vector  normal  to  surface 

tension  in  belt 

magnitude  of  force  along  n 

coefficient  of  friction  along  T 

coefficient  of  friction  along  T(g)n3  ( © ., 
vector  cross  product) 

differential  length  of  belt 

angle  between  n and  N. 


We  have  the  equation: 


' - fn  - vfT  - vf  (T®n) 


(6.9) 


From  the  properties  of  surface  curves.  Reference  4,  page  283,  we 


have 


^7  = KN,  where.’  K is  the  curvature  of  the  belt 

N = n cos  <j>  (T(g)n)  sin  <p 

T©N  = B = n .sin  <f>  + (T@n)  aos  <j> 

From  Equations  6.9,  6.10  and  6.11 


= 

ds 


-vf 


(6.10) 

(6.11) 


(6.12) 
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K sin  <f>  = of/qj  the  geodesic  curvature  of  the  belt  y (6.13) 

K aos  p = f/q3  the  normal  curvature  of  the  belt  k (6.14) 

Thus,  from  6.13  and  6.14 
tan  p - u . 

We  make  the  following  observations: 

1.  If  the  coefficient  of  friction,  u,  is  zero  along  the 
belt  line  the  tension  in  the  belt,  q3  is  a constant. 

2.  If  the  coefficient  of  friction,  u,  perpendicular  to  the 
belt  line  is  constant,  the  angle,  <j> , between  the  normals 
to  the  surface  and  the  belt  is  constant.  If  u is  zero 
the  normals  are  the  same  and  the  belt  is  a geodesic. 

Thus,  for  a frictionless  belt: 


* = o and  (6.15) 

K = f/q. 

The  equation  of  a belt  on  a nondeformable  surface  may  be  derived 
as  follows: 

Let  r = vector  from  origin  to  point  P on  the  belt.  We  have 
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dr 

Id 


kn 


y (T  © n) 


(6.16) 


= KN  = 


dn 

ds 


T-  grad  n = -k.T  + t(T  (x)  n) 


(t  is  the  geodesic  torsion) 

At  the  point  P,  r and  T are  known,  n and  grad  n are  functions  of  the 
equations  describing  the  surface.  The  normal  curvature  k may  be  computed 
from  n as: 


k = 


-T  • 


dn 
ds  ' 


From  Equations  6.13  and  6.14,  y = u k.  Hence  Equation  6.16  becomes: 


dT  d2r 

ds  j 2 


o k ( T @ n) 


(6.17) 


Equation  6.17  is  the  equation  for  the  belt.  For  a unit  sphere  we  have: 
n = -r 

grad  n = -I  (the  identity) 

t 

ds 

k-  13  normal  curvature  of  all  curves  on  a unit  sphere  is  1 

t = 03  geodesic  torsion  is  zero,  all  curves  on  a sphere  are 
lines  of  curvature. 
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Equation  6.17  for  a sphere  is: 


+ u r = 0 (6.18) 
ds 


Crossing  this  equation  with  r yields: 


— d2r  dr  dr , , 

r ® TT  u (r  ' ds  ~ ds  = 0 

ds 


but 


i2- 

— d r 
r © 2 

ds* 


d — 

^®ds^  and  r • ^7  is  zero  for  the  sphere.  Hence, 


d A dr , 


dr 


ar  fr@3a;  ~ u a? 


(6.19) 


For  constant  u.  Equation  6.19  may  be  integrated  to  give 


r © 

li 

o (r-rQ) 

+ r0  © 

r © 

li 

1 sT 
\l 

o 

+ 70®Tt 

dr. 


(6.20) 


Substituting  6.20  into  6.18  yields: 


2 2 

Letting  K = 1 + u = r n and  noting  that 


T = Tn  at  s = 0, 


solution  of  Equation  6.21  yields: 


r = ~ [u2  + aos  Ks ] + TQ  — ~ rQ  @ TQ  (aos  Ks-1) 

K K 


(6.22) 


J—  V Q 

T = -j—  = - — sin  Ks  + Tn  aos  Ks  - rn(x)  T.  sin  (Ks) 
as  K 0 K 0 ^ 0 


(6.23) 


— dT 
KN  = ^-  = 
ds 


- r g aos  Ks  - K sin  Zs  - u @ aos  Ks 


(6.24) 


From  Equation  6.14,  f = kqt  and  from  Equation  6.12 


= -\ikq  = -\xq  (k  = 1 on  unit  sphere) 


If  y is  a constant,  this  equation  may  be  integrated  to  give 
q = q0e'VS  (6.25) 


Thus,  for  the  case  of  a belt  lying  on  a unit  sphere  with  constant 
coefficients  of  friction.  Equations  6.22,  6.23  and  6.24  describe  the  geometry 
and  Equation  6.25  gives  the  tension  as  a function  of  the  arc  length  s. 


To  complete  the 
length  of  the  constrained 
be  q = F (strain) , where 
have 


analysis,  let  L(s)  be  the  variable  describing  the 

beit.  Let  the  stress-strain  relation  of  the  belt 
ds 

the  strain  is  If  the  function  is  linear,  we 


q = a 


ds 

dL 


where 


6 F 

6 strain 


is  a constant. 
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/ 


m dL  a a qa 

Thus,  tt-  = — = — 

’ ds  q qQ 


L = 


a (e^S  - 1) 


(6.26) 


or 


7 ^ n 

s = - Jln[2  + — - L] 
u a 


where 


L = 0 for  s = 0 . 


We  may  rewrite  Equation  6.26  in  terms  of  q as: 


L = — [—  - 1 L 

v q 


or 


q = qQ/(l  + y — L) 


W0 

If  the  quantity  L is  small  compared  to  one,  we  have  from 


Equation  6.26: 


s ~ — L or  qn  = as/L 
a ^ 0 
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Assume  the  belt  lies  in  a plane  and  the  surface  is  circular  (spherical) 
From  the  geometry: 


• „ 2 72  2 , 2 ,2. 

stn  d = f<p  d a + b - r + ar)/ (a  + b ) 


(6.27) 


cos  0 = (a  ^ o'  + b^  - - br)/(a‘  + b%) 

The  total  length,  Lj  of  the  belt  is: 

L = 2 [a  sec  0 + r ( 0 - tan  0)] 

= 2 [r0  + ^ + b^  - v 2 ] 

Let  q = Z (L/Lq  - D > tension  in  belt, 

/ = k (r 0 - r) j normal  force  per  unit  length 

Z = strain  coefficient,  lbs.  per  unit  stress 

k = force  coefficient  for  surface  (lbs.  per  inch  length 
per  inch  deformation) 

Tg  = undeformed  radius 

From  Equation  6.15  we  have: 

q = vfy  or 

Z (^Lq  - 1)  = kr  (v g - v) 


(6.28) 


(6.29) 
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In  the  case  where  the  surface  is  rigid  (k  = °°)  Equation  6.29  is 
not  applicable;  the  value  of  r does  not  change  and  L is  computed  from 
Equation  6.28  as  a function  of  b.  When  the  surface  is  not  rigid  (k  < °°) 
Equations  6.27,  6.28  and  6.29  must  be  solved  simultaneously  to  evaluate 
L,  0 and  r as  a function  of  b. 


Consider  the  strains 


a = L/L„  - 1 


ds_  _ 1_  6L_  dr_  6£ 

db  Lq  6r  db  6b 


but 


6L 

— = 28  and  from  Equation  6.29 


Hence, 


a ds.  = v \r  _ 2rl  — 

1 db  K lr0  db  3 


ds 

db 


SL/ 


6b 


1 + 2l8/(kLQ(2r  - rQ)) 


dq  _ „ ds  _ 


db 


* SL/Sb 


db 


1 + 218/ (kL Q (2r  - v Q) ) 


(6.30) 


Thus,  the  effective;  slope  of  the  stress-strain  curve  is  reduced  by 
the  factor  in  the  denominator  of  Equation  6.30.  Note  that  for  k = °°  (non- 
deformable  surface)  the  denominator  is  1. 
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Numerical  Example 


Consider  a case  where  a = = 7",  b = 20".  We  have  0 = it /2 

and  when  v = r 

o 

1 + 2lQ/(k  Lq  (2r  - vQ))  = 1 + ^ ? (4Q  + ?^) 

If  a belt  produces  a stress  (tension)  of  1000  lbs.  for  a strain  of 
0.1,  then  i = 10000.  If  the  7"  radius  is  reduced  to  6"  by  the  1000#,  then 
k = 1000/6  (Equation  6.29).  Thus,  i/k  = 60  and  the  slope  of  the  stress- 
strain  is  reduced  by  the  factor  1.43.  Failure  to  allow  for  the  compliance 
of  the  surface  would  produce  a 30%  error  (yjj  = i-n  t^ie  computed  stress. 

As  b increases,  r is  reduced  and  the  reduction  factor  becomes  even  greater. 

6.3.2  Air  Bag  Analyses 

(a)  Technique  for  Generating  a Set  of 
Equally  Spaced  Points  on  a Sphere 

Assume  that  along  any  two  great  circles  of  a sphere  we  are  given 
a set  of  equally  spaced  points,  X ^ 3 and  X ^ 3 as  illustrated  in  the 
sketch  below. 


X/l  *zi  *57 
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Let  the  spacing  (arc  length)  be  £ and  the  angle  between  the  vectors 
from  the  origin  to  any  two  adjacent  points  be  A = l/r}  where  r is  the  radius 
of  the  sphere. 

We  wish  to  complete  the  array  Z^  such  that  the  spacing  of  any  two 
adjacent  points  is  £.  This  can  be  done  by  starting  with  X ^ and  Z^  and 
finding  X Then  from  X ^ and  X31  finding  X^2j  etc- 


This  may  be  done  by  using  the  relation: 

X = -r~- [cos  A (X  + X ) + J( 

1 + oos  6 12  1 szn  S/2 


_ 2 

where  aos  6 = Z^  . Z^/r  . 

Letting  XJ  = X^  X 2 = X^  then  X is  Z^. 

Consider  the  segment  of  the  surface  of  the  sphere  depicted  in  the  following 
sketch  where  X+  is  the  value  of  Z using  the  plus  square  root  and  Z 

is  the  value  for  the  negative  sqpare  root. 


We  have 

cos  6 = * . X2/r2 

a v"  / 2 , 2,6,  „ 2 6, 

aos  9 = X+  . X /r  = tan  (—)  aos  2 A sea  (-^) . 


C9 


\ 
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Let  1 1 = a.X i + yX^  @ 


We  have  t ^ . X ^ = 0 = aX ^ . X ^ + QX^  . X^  = 0 


We  note  from  Equation  6.31  that 


cos  A Xj  . Xg  + sin  A . X^  © X^  = aos  A X^  . X^ 


Hence,  t . ^ @ = cos  A fZ,,  . Xp  - X^  . Xp)  = y Z? © Z,,  . Z? @ Z 


2 2 1 * 2'  ' 1 ^ 2 1 ^ 2 


Thus, 


X . X 

*i  =B  (x2-  =~ — x2)  + 

xi  ■ Ki 


aos  A fZ9  . Z9  - Z7  . Z9) 

- --J -*,©*2 


Since  ^1  ' = we  have 


8 sin  A = 


1 + 


-2 


Collecting  terms  we  obtain 


X = -t— — ^ [cos  A fl,  + XJ  + 

1 + aos  6 12 


. . . x 2 , Z9-©Z7 

fstn  A \ -1  2^  1 

sin  6/f 
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where 


X2  . X2  = X2  . X2  = r2  and  cos  6 = X2  . X2/r2 

Note  that,  for  a solution  to  exist,  sin  A sin  6/2,  i.e.,  the  arc  length  between 

the  points  must  not  be  greater  than  2 Ar. 

(b)  Effect  of  Bag  Stretch 

The  sketch  below  depicts  a half-cross  section  of  the  contact  of  a 
penetrator  of  radius  r ^ with  a bag  of  radius  r^.  The  deformed  bag  shape  is 
assumed  to  consist  of  two  circular  arcs  A-B  and  B-C.  If  the  bag  does  not 
stretch,  the  length  of  line  A-B-^C  is 
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Assume  the  bag  stretches  by  a factor  p over  the  arcs  A-B-C. 

Then, 

U2\d > = r <j)  + (r*  - p)  (<p  + \p) 

b a b 

p qos  ip  + (r  + - p)  cos  <p  = d (6.32) 

p sin  i p - (vq  + - p)  sin  <p  = 0 

The  unknowns  are  p.,  <f>  and  ip . 


If  p = 1 (no  stretch)  the  solution  reduces  to 

d 


<j>  = \pj  qos  4>  — 


r + r,  J 
a b 


P = 


v + r, 
a b . 


(r  + r, ) svn  <p  d - (r  + r,)  qos  4> 
...  , a b a b 

We  have  p = — - 


szn  \p  + szn  4> 


qos  \p  - qos  <p 


sin  ip  - sin  <p 


sin  <p 


(r  + rh)2-d2 
a b 


2 2 

d +(r  + r,)  -2d(r  + r,)  qos  <J> 

a b a b 


(r  + r,)2-d2 
a b 


2 2 

2(r  + r,)  (r  + r,-d  qos  <p)-((r  + r,)  -d  ) 

a b a b a b 


p = r + ru  - 
a u 


f , ,2  ,2 
(r  f r*)  -d 
a b 

b 2\v  t r,-d  qos  4> J 
a b 


[(v-Dr^  + p]  ip  = (?a  + - p)  <p 


These  equations  may  be  solved  using  an  iterative  type  procedure. 
As  p increases  p ->  r^.  This  is  a limiting  case  for  this  model.  At  this 
limit : 
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\ 


i 

'p 


sin 
sin  ip 


aos 


d 


2r 


a 


2r  d 
a 


Z 


In  a cylindrical  coordinate  system  (see  sketch  above)  the  surface 
of  the  bag  may  be  described  as : 


2 2 2 

r + z = r,  for  v > r,  sin  ip 
b —b 


2 2 2 

r + (z-d)  = v for  v < v sin  cb 

a —a 


2 2 2 

(r-p  sin  \pj  + (z-p  aos  i p)  = (rh-p)  3 for  sin  <p  <_  r <_  sin  ip 


where  x = v aos  03  y = r sin  9 
and  the  penetration,  p = - d 

Assume  the  bag  does  not  stretch  in  the  planes  Q = 0 and  6 = tt/2 
which  represent  the  warp  and  fill  directions  of  the  fabric.  Let  the  stretch 
factor  be  represented  by  the  first  term  of  a Fourier  Series  as 
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|i  = / /■  ( ) (I  - C()i',  (4  <>))  ~ I f |'.  II  > 


Thus,  pj  and  <\>  determined  by  Equations  6.32  are  functions  of  0. 

Sample  calculations  based  on  the  foregoing  analysis  were  performed 

to  illustrate  how  the  bag  volume  and  the  effective  contact  area  vary  with  the 

stretch  factor  8.  The  effective  area  is  the  area  obtained  by  projecting  the 

3 3 

curve  defined  by  — = 0 onto  the  XY  plane.  For  these  calculations  a 

4 inch  radius  rigid  sphere  was  assumed  to  penetrate  2.5  inches  into  two 

different  size  bags  having  radii  of  curvature  of  10.38  and  21.33  inches, 

* 

respectively.  The  results  of  the  computations  are  presented  in  Table  6-1. 
Note  that  the  volumes  and  the  effective  areas  decrease  significantly  as 
the  stretch  factor  increases.  The  current  version  of  the  CVS  program  does 
not  allow  for  this  effect  and  it  would  be  difficult  to  account  for  it  since 
it  is  a function  of  the  local  properties  of  the  bag  at  the  point  of  contact. 

Table  6-1 

EFFECT  OF  STRETCH  FACTOR  ON  AIR  BAG  VOLUME  AND  EFFECTIVE  CONTACT  AREA 


Stretch  Factor,  6 


Volume  ~ in. 


r.  = 10.38  vn. 
b 


r,  = 21.33  in. 
b 


Area  ~ in. 

r-,  = 10.38  in.  r,  = 21.33  in. 
b b 


0 

71.24 

139.92 

54.44 

100.60 

0.05 

62.84 

102.52 

52.36 

92.72 

0.10 

58.84 

88.72 

49.80 

83.88 

0.15 

56.72 

81.56 

47.56 

77.68 

0.20 

N/A 

77.32 

45.40 

71.24 

These  values  are  the  approximate  maximum  radii  of  curvature  of  the  ellipsoidal  bag 
used  for  the  static  tests  described  in  Volume  2 of  Reference  1. 
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6.4 


Table  1-1,  Question  10 


Question  10  of  Table  1-1  is: 

Is  it  necessary  to  preserve  similitude  of  certain  dynamic  degrees 
of  freedom  which  include  inertial  effects  for  deformations  associated  with 
contact? 

To  answer  this  question  consider  the  simple  dynamic  systems  illustrated 
in  Figure  6-4  below.  In  the  system  depicted  in  part  (a)  of  the  figure,  mass  m 
moving  with  velocity  v contacts  mass  m1  which  is  coupled  to  mass  m0  with 

JL  Cj 

a linear  spring.  Mass  m t is  coupled  to  a fixed  reference  with  a linear 
spring.  Masses  m ^ and  are  initially  at  rest.  Part  (b)  of  the  figure 

shows  an  equivalent  system  if  masses  m ^ and  m ^ are  negligible. 


(a) 

k - spring  constant 


(b) 

Figure  6-4  SPRING  MASS  CONTACT  MODELS 


2k 


m 


m. 


m2 

2k 


96 


Assume  that  the  initial  contact  is  sufficiently  'soft'  so  that  the 
effects  of  an  impulsive  contact  need  not  be  considered.  (Impulsive  contacts 
are  discussed  in  Section  6.2  of  this  report.)  Further,  assume  that  mass  m 
combines  with  mass  m and  that  they  move  as  a single  unit.  The  initial 

* 171  / * 

velocity  is  then  (from  conservation  of  momentum)  V-  = V m where 
* 1 
m = m + m^. 


The  equations  of  motion  are: 


m*X  = 2k  (X2  - X - (X2  - X)  ) 


m2X2  - -*  (X2  ' XV  + 2k  (X3  - X2  - <XZ  - X2)0> 


where:  X = position  of  combined  mass  m 


X^=  position  of  mass  m, 


C7s  position  of  fixed  reference. 

O 


Solving  these  equations  for,  X3  the  acceleration  of  X3  yields: 

+ V-  + p/  ^4+p2] 


k_  T* 
; 

m 


X = -^V  [[J- 

\4+ 


p sin  t 


] 


sin  u 


] (6 


2 k 

where  : = — j 

m 


2 + p + \4+  p 


2 k 2 + p 


+ V4+p‘ 


m 


p = m2/m 


.33) 
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The  solution  for  the  system  shown  in  Figure  6-4 (b)  may  be  obtained 

* 

from  Equation  6.33  by  letting  p = 0 and  m = m.  Thus  for  this  system: 


k_  y sin  u )t 
ml  a) 


V \ I — sin  (t  tf  —J 

1 y m y m 


(6.34) 


where  co  = k/m. 

We  only  show  the  solution  for  the  acceleration  of  X since  in  the 
crash  environment  most  injury  criteria  are  based  on  accelerations  or  forces. 
For  simplicity  let  us  normalize  Equation  6.33  by  dividing  by  the  amplitude 
factor  of  Equation  6.34.  We  have 


a = X/(-~  VjuJ 
mi  l 

* 

= t[  Ul  - p/V 4+p2> 

= A - sin  + A0  sin 


With  the  same  normalization  Equation  6.34  becomes: 

a.Q  = sin  ut  (6.36) 


— s^n  to 

Uj  1 


t + (1  + p/  y4+ si 


svn  to ^t\ 


(6.35) 


For  this  simple  model  Question  10  may  be  rephrased  as:  Is  the  difference 

between  Equation  6.35  and  Equation  6.36  significant?  We  cannot  give  an  absolute 
answer  to  this  question  since  the  difference  of  results  depends  on  the  problem 
at  hand.  We  do,  however,  present  numerical  results  showing  how  the  differences 
depend  on  the  masses  and  spring  constant. 

The  numerical  results  are  given  in  the  form  of  tables.  The  first  table, 
Table  6-2,  shows  how  the  frequency  and  period  of  a simple  spring  mass  system 
(Equation  6.34)  varies  with  the  ratio  of  the  spring  constant  to  the  weight.  The 
formulas  are: 
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Table  6-2 

FREQUENCY  AND  PERIOD  OF  SIMPLE  SPRING  MASS  SYSTEM 


Spring  Constant/Weight 

K/W  (in.  -1) Freq  (HZ) 


Period  (MS) 


1. 00D-01 

0.99 

1011.31 

1 . 78D-01 

1.32 

758.38 

3. 16D-01 

1 . 76 

568.70 

5.62D-0I 

2.34 

426.47 

1 - 00D+00 

3. 13 

319.81 

1 . 78D+00 

4.  17 

239.82 

3. 1 6D-H*  0 

5.56 

179.84 

5.62D+00 

7.42 

134.86 

1.00D+01 

9.99 

101.13 

1.78D+01 

13.19 

75.84 

3. 16D+01 

17.58 

56  ■ 8 r' 

5.62D+01 

23.45 

42. 65 

1. 00D-KI2 

31.27 

31.98 

1 - 78D+02 

41.70 

23.98 

3.  16D-KI2 

55.60 

17.98 

5.62D+02 

74.  15 

13.49 

1 - OOD+-03 

98.88 

10.11 

1 . 78D-K»3 

131.86 

7.58 

3. 16D+03 

175.84 

5.69 

5.62D+03 

234.48 

4.26 

1. 000+04 

312.69 

3.20 

1 .78D+04 

416.98 

2.40 

3. 160+04 

556. 05 

1.80 

5. 620+04- 

741.50 

1.35 

1. 000+05 

988.81 

1 . 01 
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frequency  f 


1_ 
2 TT 


0) 


1_ 

2n 


9 , Hertz 


period  T 

where  W 

k 

9 


1000/ 


f j milliseconds 


mg y pounds 


spring  constant  pounds /inch 


acceleration  of  gravity 3 386  inches /second11 


From  Table  6-2  we  note  that  for  k/w  = 10  the  period  is  about  100 
milliseconds.  The  ratio  of  k/w  must  be  100000  to  give  a period  of  1 milli- 
second. 


The  second  table.  Table  6-3,  shows  the  effect  when  mass  2 of  the 
3-mass  contact  model  is  negligible  (p  = 0 ) but  mass  1 is  not.  The  formulas 


■ 

* t 

= I''*  + velocity  ratio 


Aj  = 1/(1  + r)^^y  amplitude 


Tj/T  = J 1 + r y period  ratio 


r = m^/m  , mass  ratio 

JL 

The  third  table.  Table  6-4,  shows  the  response  of  the  3-mass  system  when  mass  1 

is  negligible  and  p (the  ratio  of  to  m ) is  varied.  The  formulas  are: 

2 
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Table  6-3 


EFFECT 

OF  MASS  RATIO  ON  RESPONSE 

OF  3 -MASS  CONTACT  MODEL 

WITH  NEGLIGIBLE 

MASS  2 

Mass  Ratio 

Velocity  Ratio 

Period  Ratio 

Amplitude 

Mj/M 

Vl*/V 1 

t1  n 

Al 

1 . 000+01 

0.  0909- 

3.31 66 

0. 0274 

5.620+00 

0. 1510 

2.5736 

0. 0587 

3. 160+00 

0.24-»)3 

2 . 0402 

0. 1178 

1 . 78D+do 

0. 3599 

1.6668 

0.2159 

1 . 000+00 

0.5000 

1.4-142 

0.3536 

5.62D-01 

0.64-01 

1 . 2499 

0.5121 

3. 160-01 

0.7597  ! 

1.1473 

0.6622 

1 .78D-01 

0.84-90 

1 . 0853 

0.7823 

1 . 000-01 

0.9091 

1 . 0488 

0. 8668 

5.62D-02 

0.9468 

1 . 0277 

0.9212 

3. 160-02 

0.9693 

1. 0157 

0.  9544 

1.780-02 

0.9825 

1 . 0089 

0.9739 

1 . 000-02 

0.9901 

1 . 0050 

0.9852 

5.62D-03 

0.9944 

1 . 0 028 

0.9916 

3. 160-03 

0.  9968 

1 . 0016 

0.9953 

1.780-03 

0. 9982 

1 . 0009 

0. 9973 

1 . 000-03 

0.9990 

1 . 0 0 05 

0.9985 

5 . 620— 04> 

0.9994 

1 . 0003 

0 . 9992 

3.  160—04 

0.9997 

1.0002 

0.9995 

1 . 780—04- 

0.9998 

1 . 0001 

0. 9997 

1 . 000-04- 

0 . 9999 

1 . 0 0 0 0 

0 . 9999 

5.620-05 

0 . 9999 

1.00 0 0 

o . 9949 

3. 160-05 

1.00  0 0 

1 . b 0 0 0 

1 . 0000 

1 . 780-05 

1 . 0 0 0 0 

1.00 0 0 

1 . 0000 

1 . 000-05 

1 . 0000 

1. 0000 

1 . 0 0 0 0 
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Table  6-4 


EFFECT  OF  MASS  RATIO  ON  RESPONSE  OF  3-MASS  CONTACT  MODEL 
WITH  NEGLIGIBLE  MASS  1 


Mass  Ratio 


Period  Ratios 


Amplitudes 


Amplitude 

Ratio 


Mr/M 

T2/T 

t2/t 

VT2 

Ai 

- A2 

a2/ai 

1 . 00D+01 

2.3557 

0.6712 

0.2849 

0. 0457 

1 . 3293 

29. 0585 

5.62D+00 

1 . 8434 

0.6432 

0.3489 

0 . 1 066 

1 . 2493 

11.7219 

3. 16D+00 

1.4820 

0.5959 

0. 3994 

0.2310 

1 . 0996 

4. 7597 

1.78D-KI0 

1 . 2703 

0.5249 

0.4132 

0.4262 

0. 8737 

2.  0498 

1 . 00D+00 

1.1441 

0.4370 

0.3820 

0 . 6325 

0.6325 

1 . 0000 

5. 62D-01 

1 . 077  0 

0.3481 

0.3232 

0.7855 

0.4424 

0. 5632 

3. 16D-01 

1.0418 

0.2699 

0.2591 

0.8791 

0.3121 

0.3550 

1 . 78D-01 

1.0230 

0.2061 

0.2015 

0.  9324 

0.2244 

0. 2407 

1 . 00D-01 

1.0127 

0. 1561 

0. 1542 

0.9622 

0. 1639 

0.1704 

5.62D-02 

1 . 0 07 1 

0.1177 

0. 1169 

0. 9788 

0. 121 0 

0. 1237 

3. 16D-02 

1. 0040 

0. 0886 

0. 0882 

0.9881 

0. 0900 

0. 0910 

1 . 78D-02 

1 . 0 022 

0. 0665 

0 . 0664 

0. 9933 

0. 0671 

0. 0676 

1 . 000-02 

1. 0013 

0. 0499 

0.0499 

0 . 9962 

0 . 05  02 

0. 0504 

5.620-03 

1.0007 

0. 0375 

0. 0374 

0. 9979 

0. 0376 

0. 0377 

3. 16D-03 

1 . 0004 

0. 0281 

0. 0281 

0 . 9988 

0.0282 

0. 0282 

1.78D-03 

1 . 0002 

0. 0211 

0. 0211 

0 . 9993 

0. 0211 

0. 0211 

1 . 0 OD— 03 

1 . 0001 

0. 0158 

0. 0158 

0 . 9996 

0. 0158 

0. 0158 

5. 620—04- 

1 . 0001 

0.  0119 

0. 0119 

0 . 9998 

0.  01 19 

0. 0119 

3. 16D-04 

1. 0000 

0. 0089 

0. 0089 

Q _ QQQQ 

0. 0089 

0. 0089 

1.78D-04 

1 . 0000 

0. 0067 

0. 0067 

0 . 9999 

0. 0067 

0 . 0 067 

1 . 00D— 04- 

1 . 0000 

0 . 0 05  0 

0. 0050 

1 . 0000 

0 . 0 05  0 

0. 0050 

5.62D-05 

1.0000 

0. 0037 

0. 0037 

1 . 0000 

0. 0037 

0. 0037 

3. 16D-05 

1 . 0000 

0. 0028 

0. 0028 

1. 0000 

0. 0028 

0. 0028 

1 . 780—05 

1 . 0000 

0. 0021 

0. 0021 

1 . 0000 

0. 0021 

0. 0021 

1 . 00D-G5 

1.0000 

0. 0016 

0 ~ 0 0 1 6 

1.0000 

0. 0016 

0. 0016 
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T/T  = 


vr  = 


y (2  + 0 + J 4+p2)/4  , 

y 0/(2  + P + ^ «+p2;  , 


(1  - p/ y4+p  ) (o/w^  j 

fi  + p/  ^ 4+p2)  u)/(jj2  , 


period  ratio 


period  ratio 


amplitude 


amplitude 


The  ratios  of  T2//Tl  anc*  ^2^1  are  a*so  tabulated. 

In  conclusion  the  answer  to  the  question  as  to  whether  or  not  it  is 
necessary  to  include  inertial  effects  for  contacts  is  dependent  on  the  situation. 
In  this  section  we  give  numerical  results  for  a simple  contact  model  which 
should  be  helpful  to  the  program  user  in  making  a decision  as  to  whether  or 
not  inertial  effects  should  be  included.  Table  6-2  which  shows  the  periods  of 
oscillation  as  a function  of  stiffness  to  weight  is  primarily  of  use  in 
estimating  the  integrator  step  size  or  in  reaching  a decision  on  the  use  of  an 
impulsive  contact.  The  numerical  integrator  should  have  a maximum  step  size 
that  is  no  greater  than  1/10  of  the  period  (this  is  relative  to  the  accuracy 
desired);  thus.  Table  6-2  shows  that  if  one  wishes  to  use  a maximum  step  size 
of  1 millisecond  the  stiffness  to  weight  ratio  should  be  no  greater  than  1000. 
Table  6-3  shows  the  effects  of  ignoring  part  of  the  system  mass  such  as,  for 
example,  the  mass  of  the  material  which  is  carried  with  the  head  on  a windshield 
contact.  For  a mass  ratio  of  1/10  there  is  a 10%  change  in  velocity,  a 5% 
change  in  period  and  a 13%  change  in  amplitude.  Ta^le  6-4  shows  the  effects 
of  the  first  two  modes  of  vibration  of  an  elastic  impact.  In  an  actual  physical 
problem,  since  we  are  dealing  with  continuous  structures,  there  will  be  an 
infinite  number  of  modes  of  vibration.  In  an  impact  problem  the  higher  order 
modes  will  be  of  progressively  less  importance;  that  is,  these  amplitudes 
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will  decrease  as  the  frequency  increases.  Table  6-4  shows  that  if  the 
effective  mass  of  the  impacted  structure  (M^)  is  1/10  of  the  effective 
impacting  mass  {M)  the  period  of  the  principal  resonance  Tj  differs  by 
about  1%  from  the  period  that  would  occur  if  mass  2 were  ignored.  The 
secondary  resonance  has  a period  that  is  16%  of  the  primary.  The  amplitude 
of  the  principal  resonance  differs  by  4%  and  the  secondary  resonance  has  an 
amplitude  of  17%  of  the  primary.  The  principal  effect  of  ignoring  mass  2 
would  be  the  neglection  of  the  17%  amplitude  secondary  frequency.  It  should 
be  remembered  that  these  values  pertain  to  the  acceleration  response;  the 
effects  on  position  will  be  reduced  by  the  square  of  the  period  ratio  (T^/T^) 


i 
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7.0 


CONCLUDING  REMARKS 


Although  this  report  may  seem  to  the  reader  to  be  a collection  of 
unrelated  topics,  the  topics  have  one  thing  in  common,  i.e.,  they  are  all 
related  to  various  aspects  of  the  CVS  program. 

The  user  must  remember  that  the  CVS  program  is  based  on  rigid  body 
dynamics  and,  hence,  when  used  to  simulate  real-life  conditions  can  only 
approximate  the  true  dynamics  of  the  system.  When  the  program  is  used  to 
simulate  a dummy,  the  skeletal  structure  of  the  dummy  can  be  modelled  quite 
accurately  but  the  presence  of  the  'skin,'  and  the  rubber  neck  and  spine 
preclude  the  precise  modelling  with  a program  based  on  rigid  body  dynamics. 

The  work  on  the  shoulder  model  reported  in  Section  3.2  of  this  report 
demonstrates  the  flexibility  of  the  program  in  allowing  the  user  to  define  a 
more  precise  model  of  a dummy.  The  various  topics  discussed  in  Section  6 of 
this  report  point  out  some  of  the  limitations  of  the  model  and  attempt  to  give 
some  insight  to  the  effects  of  using  impulsive  forces  instead  of  force- 
deflection  relationships  for  "hard”  type  contacts,  the  use  of  point  contact 
algorithms  instead  of  deformable  contact  algorithms,  and  the  use  of 
simplified  models  of  belt  and  air  bag  restraint  systems  instead  of  more 
precise  models  of  these  restraint  systems. 

The  Response  Measure  Approximating  Function  Generator  described  in 
Section  4 should  prove  to  be  a useful  tool  for  the  user  who  wishes  to  perform 
parametric  studies  and  have  £i  means  of  interpolating  the  response  measures  as 
a function  of  several  parameters. 

The  dynamically  equivalent  system  algorithm  described  in  Section  5 
of  this  report  has  no  direct  bearing  on  the  CVS  program  but  points  out  a little 
known  fact  about  a system  of  interconnected  rigid  bodies.  In  particular,  it 
shows  that  the  mass  distribution  to  produce  a particular  system  response  is  not 
unique.  This  implies  that  the  exact  internal  structure  of  an  interconnected 
rigid  body  system  cannot  be  determined  from  the  observation  of  external  dynamic 
and  kinematic  responses.  The  theorem  may  prove  useful  in  defining  "canonical 
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models'  of  systems,  i.e.,  a "canonical  model'  may  be  defined  as  one  where  all 
of  the  segments  have  the  same  mass. 

Finally,  the  CVS  program  is  a useful  tool,  but  the  success  of  its 
application  to  a particular  problem  requires  a full  understanding  of  its 
capabilities  and  limitations  by  the  user.  There  is  much  work  yet  that  could 
be  done  to  improve  the  program,  particularly  the  development  of  new  contact 
algorithms  to  better  model  the  interaction  of  deformable  bodies  which  is  the 
real  world  situation. 
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APPENDIX  A 


PLOTS  OF  RESPONSE  MEASURES  FROM  SHOULDER  MODEL  SIMULATION  STUDY 


A-l 


This  appendix  contains  plots  of  selected  response  measures  for 
the  six  CVS  runs  performed  as  part  of  the  study  of  dummy  shoulder  models 
that  are  described  in  Section  3.2  of  this  report.  For  each  run  the  plots 
are  presented  in  the  following  order: 

(a)  Right  lower  arm  resultant  acceleration 
versus  time 

(b)  Right  lower  arm  Z displacement  versus  X displacement 

(c)  Right  lower  arm  Y displacement  versus  X displacement 

(d)  Shoulder  yoke-clavicle  pivot  flexure  angle  versus  time 

(e)  Precession  and  nutation  angles  of  the  shoulder  Euler 
joint  versus  time 

(f)  Precession  and  nutation  angles  of  the  elbow  Euler 
joint  versus  time 

(g)  Upper  torso  resultant  acceleration  versus  time 

(h)  Upper  torso  pitch  angle  versus  time. 
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Figure  A-l  RUN  NO.  1 RESPONSE  MEASURE  PLOTS 
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A-6 


W RY  FLEXURE  VS.  TIME 
UNLOCK.  P _ LOCK.  N - UNLOCK.  S r LOCK 

Figure  A-l  (Cont'd.) 


A- 7 


C6 


09 


— 1 o 
06- 


(e)  RS  PRECESSION  AND  NUTATION  VS.  TIME 
UNLOCK,  P = LOCK,  N - UNLOCK,  S - LOCK 
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(f)  RE  PRECESSION  HND  NUTHT I ON  VS.  TIME 
UNLOCK.  P - LOCK.  N = UNLOCK.  S = LOCK 
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(g)  UPPER  TORSO  RESULTANT  ACCELERATION  V S.  TIME 
: - UNLOCK  , P - LOCK,  N = UNLOCK,  S _ LOCK 

Figure  A-l  (Cont'd.) 
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(h)  UPPER  TORSO  PITCH  VS.  TIME 
UNLOCK.  P = LOCK.  N = UNLOCK.  S = LOCK 
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UNLOCKED. r _ LOCKED. P.N  r UNLOCKED. S - LOCKED 

Figure  A-2  (Cont'd.) 
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(e)  RS  PRECESSION  AND  NUTATION  VS.  TIME 
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Figure  A-2  (Cont'd.) 
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(g)  UPPER  TORSO  RESULTANT  ACCELERATION  VS-  TIME 
X.Z  = UNLOCKED , T = LOCKED. P.N  - UNLOCKED. S - LOCKED 
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(h)  UPPER  TORSO  PITCH  VS.  TIME 
UNLOCKED. Y = LOCKED. P.N  = UNLOCKED. S = LOCKED 

Figure  A-2  (Cont'd.) 
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Figure  A- 3 RUN  NO.  3 RESPONSE  MEASURE  PLOTS 
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(d)  RY  FLEXURE  VS.  TIME 
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(e)  RS  PRECESSION  AND  NUTATION  VS  . T INF- 
UNLOCKED  . Z * Y = LOCKED.  P.N  = UNLOCKED.  S - LOCKED 


II 

X 


A-24 
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(h)  UPPER  TORSO  PITCH  VS*  TIME 
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Figure  A-4  (Cont'd.) 
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Figure  A- 5 RUN  NO.  5 RESPONSE  MEASURE  PLOTS 
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(f)  RE  PRECESSION  RND  NU  TRT I ON  VS.  TIME- 
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-D 


3 

CM 


O 

o 


o 

O' 


o 

00 


o 

CO 


o 


o 

<*> 


o 

C4 


o 


cn 


LU 


CO 

> 

2: 

CO 

a 

UJ 

t— 

iC 

i— 

C_> 

a: 

a 

\ 

oc 

_i 

• 

UJ 

X) 

— 1 

ii 

+-* 

uJ 

c 

C_) 

CO 

0 

co 

u 

CE 

* 

CO 

h- 

UJ 

LO 

Z 

1 

0 

1 

0 

0) 

1 

_i 

u 

co 

2: 

3 

CO 

CO 

W) 

•H 

uj 

U. 

QE 

! 1 

CO 

2! 

CO 

• 

0_ 

5 

r— 

• 

CO 

c ✓ 

uj 

UJ 

Q_ 

0 

CL. 

0 

CO 

— 1 

, „ 

it 

00 

v — / 

>- 

rvi 

X 


A-41 


j 


o 

U) 


o CO 
oo  s: 


UJ 
o z: 


o 

iO 


o 

C4 


tl 

>- 

rvj 

X 


A-42 


M UPPER  TORSO  PITCH  VS.  TIME 
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Figure  A- 6 RUN  NO.  6 RESPONSE  MEASURE  PLOTS 
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RY  PLFXURC  VS.  Tint: 

LOCK.  Y - UNLOCK.  P = LOCK.  N = UNLOCK.  S r LOCK 
Figure  A-6  (Cont'd. ) 
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Figure  A- 6 (Cont'd. 
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(f)  RE  PRECESSION  AND  NUTATION  VS.  TINE 
LOCK.  Y = UNLOCK.  P = LOCK.  N - UNLOCK.  S - LOCK 

Figure  A-6  (Cont'd.) 
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Cg)  UPPER  TORSO  RCSULTRNT  flCCE.l  P.RHT  I ON  VS-  TINf 
: LOCK.  Y - UNLOCK.  P = LOCK.  N = UNLOCK.  S _ LOCK 

Figure  A-6  (Cont'd.) 
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(h)  UPPER  TORSO  PITCH  VS.  TIME 
LOCK.  T - UNLOCK.  P - LOCK.  N = UNLOCK.  S - LOCK 

Figure  A-6  (Cont'd.) 


APPENDIX  B 


LISTING  OF  THE  FORTRAN  PROGRAM  FOR  THE 
RESPONSE  MEASURE  APPROXIMATING  FUNCTION  GENERATOR 


B-l 


C RMAFG  LINK  WITH  SETZ,  SETU,  SIGM,  SOLVE 

C THIS  VERSION  OF  THE  RMAFG  IS  FOR  TERMINAL  OPERATION 
C CHANGE  READ  AND  WRITE  STATEMENTS  TO  CONVERT  TO  BATCH  OPERATION 

C INPUT  (CARDS  REFER  TO  LINES  OF  DATA  SET) 

FORMAT (20A4 ) INPUT  DESCRIPTION 
FORMAT (20A4)  INPUT  DESCRIPTION .CONTINUED 


FORMAT (2014) 
NY,N,NDG(J) ,J=1,N 


NY 

N 

IF  N=0 
NDG(J) 


FORMAT (2014) 
NDATA .MODE , I VAL 


NDATA 

MODE 


IVAL 


CARD 
CARD 
CARD 
C 
C 
C 

c 
c 
c 

C 

c 

CARD 
C 
C 
C 
C 
C 
C 
C 
C 

CARDS 
C 
C 
C 

c 

CARD  NDATA+3 
C IF  IVAL 

C IF  IVAL 

C NV 

C 

CARDS  NDATA+4-NDATA+3+NV 

C 
C 
C 

c 
c 
c 
c 

£»»*»*• 

c 

c . 
c 
c 

£•*«««* 

c 
c 
c 
c 
c 
c 
c 
c 


NUMBER  OF  DEPENDENT  VARIABLES 
NUMBER  OF  INDEPENDENT  VARIABLES 
PROGRAM  TERMINATES 
MAXIMUM  DEGREE  OF  VARIABLE  J 
NDG(J-1 ) .LE.NDG(J) 


NUMBER  OF  SETS  OF  DATA  POINTS 
0 NORMAL 

REMOVE  MEAN  FROM  INDEPENDENT  VARIABLES 
REMOVE  MEAN  AND  DIVIDE  BY  STANDARD  DEVIATION 
COMPUTE  FIT  ONLY  AT  POINT  X(  ) USED  TO  EVALUATE  C 
READ  ADDITIONAL  POINTS  X(  ) TO  COMPUTE  FIT  AND  EVALUATE 


5 -NDATA +2  FORMAT ( 10F8. 0) 

(X(J,K) , J = 1 ,N) , (Y(K,L) ,L=1,NY) 

X(J,K),J=1,N  VALUES  OF  INDEPENDENT  VARIABLES  AT  K 
Y (K.L) ,L= 1 ,NY  VALUES  OF  DEPENDENT  VARIABLES  (FUNCTION)  AT  K 


FORMAT (2014) 

= 0 THIS  IS  THE  SAME  AS  CARD  1 
< 0 OR  > 0 PROGRAM  READS  VALUE  OF  NV 
NUMBER  OF  POINTS  TO  BE  EVALUATED  USING  COMPUTED  FIT 


FORMAT ( 10F8. 0) 
(X(J,M),J=1,N),(Y(M,L),L=1,NY) 
X(J,M) , J=1 ,N 
Y(M,L) ,L=1 ,NY 


ADDITIONAL  VALUES  OF  INDEPENDENT  VARIABLES 
VALUES  OF  DEPENDENT  VARIABLES  (IF  KNOWN)  AT  M 
PROGRAM  WILL  COMPUTE  YP(M)  FOR  EACH  L 


PROGRAM  STOPS 


PROGRAM  COMPUTES  THE  BEST  LEAST  SQUARE  FIT  FOR  EACH  DEGREE 
FROM  0 TO  MAXIMUM  AS  SPECIFIED  ON  CARD  3 


OUTPUT 

KM,ERR,(Z(L),C(L,I+1),L=1,KM) 

KM  NUMBER  OF  TEEMS (COEFFICIENTS)  IN  FIT 

ERR  RMS  ERROR  OF  »FIT 

Z (L ) SUBSCRIPT  OF  L'TH  COEFFICIENT 

C(L , 1+1 )VALUE  OF  COEFFICIENT  Z(L) 


B-2 


o o o o o 


PROGRAM  PROVIDES  A PRINTER  PLOT  OF  THE  HISTOGRAM  OF 
ERRORS  AND  THE  FIRST  TEN  POINTS  WITHIN  THE  SPECIFIED  ERROR  RANGE 

PROGRAM  PRINTS  THE  FUNCTION  VALUE  COMPUTED  FROM  THE  FIT 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*4  STAR, CNT,B,CNTM, CCS, BLANK, HY 
DIMENSION  HY(20,2) 

C DIMENSION  PROGRAM  FOR  3 INDEPENDENT  VARIABLES,  7 DEPENDENT  VARIABLES 
C AND  90  DATA  SETS 

DIMENSION  Y(90,7),X(3,90),YP(90),XX(3,90),XM(3),XS(3) 

DIMENSION  NDG(3).NIW(3),LIW(3) 

C DIMENSION  PROGRAM  FOR  A MAXIMUM  OF  20  COEFFICIENTS 
DIMENSION  U ( 20 ) , Z ( 20 ) , C ( 20 , 22 ) 

C DIMENSION  PROGRAM  FOR  A MAXIMUM  OF  10  POINTS  ON  HISTOGRAM 


DIMENSION 

CNT (11), MM (11), JH (11) , JHP( 10,11) 

DATA 

NMAX/3/.MMAX/  90/,KMAX/20/ 

r 

DATA 

STAR/'*'/, BLANK/'  ’/ 

w 

C 

NY 

NUMBER  OF  DEPENDENT  VARIABLES 

c 

N 

NUMBER  OF  INDEPENDENT  VARIABLES 

c 

MDEG 

MAXIMUM  DEGREE 

c 

NDG 

MAXIMUM  DEGREE  OF  INDEPENDENT  VARIABLE  J 

c 

NMAX 

STORAGE  LIMIT  ON  INDEPENDENT  VARIABLES 

c 

MMAX 

STORAGE  LIMIT  ON  DATA  POINTS 

c 

c 

KM  AX 

STORAGE  LIMIT  ON  COEFFICIENTS 

c 

Y(K,L) 

VALUE  OF  FUNCTION  L AT  POINT  K 

c 

X(J,K) 

VARIABLE  J POINT  K 

C OPEN  INPUT  FILE 

CALL  OPEN (9, ’RMAFG.DAT  ') 

2 READ(9,4)HY 
4 FORMAT (20A4) 

READ(9,6,END=12)NY,N, (NDG(J) , J = 1 ,N) 

6 FORMAT (2014) 

IF(N.EQ.0)ST0P  100 
MDEG=NDG  (N ) 

CALL  SETZ ( N , MDEG , NDG , NIW , LIW , NM , Z , NMAX , KM AX ) 

READ ( 9 , 6 ) NDATA , MODE , I VAL 

WRITE (1 , 6)NY, N, {NDG (J),J=1,N), NDATA, MODE, IVAL 
IF (NDATA. LE.MMAX)GO  TO  10 
WRITE (1 , 8)NDATA,MMAX 

8 FORMAT ('  THE  NUMBER  OF  DATA  POINTS ', 15,/ 'EXCEEDS  THE 

X ALLOWED  STORAGE' ,15,  ' PROGRAM  TERMINATED') 

STOP 

10  IF (NDATA. GE.NM(MDEG+I ))GO  TO  14 
WRITE ( 1 , 22) NDATA, NM (MDEG* 1 ) 

12  STOP 

14  DO  16  K=1, NDATA 

READ(9, 18, END=12) (XX ( J ,K) , J=1 ,N) , (Y(K, J) , J=1 ,NY) 

16  WRITE ( 1,20) (XX(J,K), 1=1, N) , (Y(K, J) , J=1 ,NY) 

18  FORMAT ( 10F8. 0)  i 

20  FORMAT (1X,F7.2,9F8.3) 

22  FORMAT ('  THE  NUMBER  OF  DATA  POINTS ',15/'  IS  LESS  THAN  THE 
X NUMBER  OF  COEFFICIENTS ' ,15,  ’ PROGRAM  TERMINATED') 
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24  WRITE (1,26) 

26  FORMAT ('  ENTER  LP,L1,L2'/ 

X ' LP  : 2 FOR  LINE  PRINTER,1/ 

X • L1,L2  SELECTS  DEPENDENT  VARIABLES,  L = L1.L2  ’) 

READ( 1 , 6)LP,L1 ,L2 
IF (LP.NE.2)LP  = 1 
DO  92  L = L 1 ,L2 

DO  28  K = 1.NDATA 

DO  28  J = 1,N 

28  X(J,K)  = XX(J,K) 

CALL  SIGM(N, NDATA ,NMAX, MODE ,YM , YS,XM,XS,Y(1,L),X) 

M=KMAX+2 
DO  30  K= 1 , KMAX 
DO  30  Js 1 ,M 
30  C(K,J)sO. 

FER=0. 

I =NM(MDEG+1 ) 

DO  36  ND=1, NDATA 
WRITE (1,32)ND 

CALL  SETU(N,MDEG,NIW,LIW,NDG,U,X( 1 ,ND) ) 

32  FORMAT ('  ND  »,I4) 

DO.  34  J= 1,1 

C(j,I+1)=C(J,I+1)+U (J)*Y(ND,L) 

DO  34  KsJ.I 

34  C(J,K)=C(J,K)+U(J)*U(K) 

36  FER=FER+Y (ND,L)**2 
DO  38  Jr  1 , 1 
38  C ( J , I +2 ) =C  ( J , I + 1 ) 

WRITE (LP,40)HY 
40  FORMAT ( 1X.20A4/20A4/) 

WRITE (LP, 42 )YM,YS 

42  FORMAT ( ' Mean  ,,F20.5,,  Sigma  ',F20.5,'  for  Function  Y'// 

X ' J Degree  ' ,5X,'  Mean',5X,'  Sigma:  for  Parameters  X(  )’) 

WRITE (LP, 44 )(J,NDG(J),XM(J) ,XS(J) , J=1 ,N) 

44  FORMAT ( IX, 13, 15, 2F 15. 5) 

IF(MODE.EQ. 1 )WRITE(LP,46) 

IF (MODE. EQ. 2) WRITE (LP, 48) 

46  FORMAT(/'  MEAN  removed  from  Parameters  X(  ).') 

48  FORMAT(/'  MEAN  removed  from  Parameters  X(  ),  X normalized.') 

12=1+2 
ArNDATA 
Mr  1 

DO  54  J=1,I 

CALL  SOLVE ( J , C , I , KMAX ) 

IF( J.LT.NM(M) )GO  TO  54 
KM=NM(M) 

M=M  + 1 

ERRrFER 

DO  50  Krl.KM 

50  ERR  rERR-C  (K,  1+1  )*C  (K,  I-r2) 

IF (ERR .GT . 0. )ERR=DSQRT (ERR/A ) 

WR ITE (LP, 52)KM, ERR, (Z(K),C(K, 1+1), K=1,KM) 

52  FORMAT(/'  The  error  using', 15,'  coefficients  is',F15.5/ 

X ' Coefficients’ , 10X, ' Value ' //25 ( IX, K 10. 0.F20.6/) ) 

54  CONTINUE 

DO  56  J = 1, 1 1 
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JH ( J ) =0 
56  CNT(J)=0. 

CNTM  =0. 

DO  60  ND=1,NDATA 

CALL  SETU (N ,MDEG , NIW , LIW , NDG , U , X(  1 , MD ) ) 

FCN=0. 

DO  58  K=1 , 1 

58  FCN=FCN+U(K)»C(K,I+1) 

C SET  UP  HISTOGRAM 

J=2.*(FCN-Y(ND,L))/ERR+6.5 

YP (ND) =FCN 

IF( J.LT. 1 )J=1 

IFCJ.GT. 11)J=11 

CNT( J )=CNT (J )+1 . 

IF(JH(J).GE. 10)GO  TO  60 
JH  (J )=JH  ( J ) + 1 
M=JH ( J ) 

JHP(M,J)=ND 

60  IF (CNT( J ) .GT. CNTM) CNTM =CNT ( J ) 

WRITE(LP, 62) 

62  FORMAT ('  Error  Distribution  -Histogram’/ 

X ' Sigma  # Points' ,20x, 'Point  IndeX') 

CCS=1. 

IF ( CNTM . GT . 20 )CCS  =20 . /CNTM 
B=-2.5 

DO  64  4=1,11 
M=CNT  (J)*CCS 

IF(M.EQ.O. AND.CNT (J).NE.0.)M=1 
M2=JH ( J ) 

Ml =21-M 

IF(M.EQ.0)WRITE(LP,66)B,CNT (J) 

IF(M.GT. 0. AND. M2. E0.0) WRITE (LP, 66)B,CNT (J) , (STAR ,K=1 ,M) 

IF(M.GT. 0. AND.M2.GT. 0)WRITE (LP , 66 )B,CNT (J),(STAR,K=1,M), 

X (BLANK, K=1,M1),(JHP(Kf J) ,K=1,M2) 

64  B=B+.5 

66  FORMAT (1X,F6. 2, F 10.0, IX, 21 A 1 , 1014) 

C PRINT  FUNCTION  Y AT  VALUES  OF  X(  ) USED  TO  DETERMINE  FIT 
WRITE (LP, 68 )L 

68  FORMAT('1  Function' ,12, ' Evaluations  at  Given  Input  Points,'/ 

X ' Point  Y(Input)  Y(Fit)  Parameters  X(  )') 

DO  72  M=1 ,NDATA 
IF(MODE.EQ.O)GO  TO  72 
DO  70  K=1,N 

IF(MODE.EQ.2)X(K,M)=XS(K)*X(K,M) 

70  X(K,M)=X(K,M) +XM (K ) 

72  WRITE (LP, 74 )M,Y(M,L) ,YP(M),(X(J,M),J=1 ,N) 

74  FORMAT ( IX, 13, 5F 12. 3) 

IF(IVAL.EQ.O)GO  TO  92 
C 

COMPUTE  AND  PRINT  EVALUATED  VALUES  AT  SPECIFIED  INPUT  POINTS 
C 

READ(9 , 6)NV 
IF (NV.EQ.O)GO  TO  92 
WRITE (LP, 76 )L 

76  F0RMAT(/'  Function ', 12, ' Evaluations  at  Specified  Input  Points.' 
X /'  Point  Y(Input)  Y(Fit)  Parameters  X(  )') 
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DO  78  M = 1 , NV 

78  READ (9 , 80 , END =92 )(X(J,M),J=1,N),(Y(M,K)fK=1 , NY ) 

80  FORMAT ( 1 0F  8.0) 

DO  90  M = 1 , NV 
IF(MODE.EQ.O)GO  TO  84 
DO  82  K=  1 , N 
X(K,M)=X(K,M)-XM(K) 

82  IF ( MODE . EQ . 2 ) X ( K , M ) =X (K , M ) /XS ( K ) 

84  CALL  SETU(N,MDEG,NIW,LIW,NDG,U,X( 1 ,M)) 

YP (M) =0. 

DO  86  K= 1 , I 

86  YP (M)  = YP(M)+U (K)*C(K, 1+1 ) 

IF(MODE.EQ.O)GO  TO  90 
DO  88  K=1 , N 

IF(M0DE.E0.2)X(K,M)=X(K,M)»XS(K) 

88  X(K,M)=X(K,M) +XM (K ) 

90  WRITE (LP, 74 )M, Y(M tL), YP (M) ,(X(J,M) , J=1,N) 

C RETURN  FOR  OTHER  DATA  SETS 
92  IF (L.LT.L2)WRITE(LP, 94 ) 

94  FORMAT (1H1) 

STOP 

END 
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COMPUTATION  OF  SUBSCRIPTS  Z 

C INITIAL  COEFFICIENT  CHECK;  MODIFIED  DEC  3 1981  FOR  WP 

SUBROUTINE  SETZ ( N ,MDEG , NDG , NIW , LIW , NM , Z , NMAX , KMAX ) 

IMPLICIT  REAL*8 (A-H ,0-Z) 

DIMENSION  NIW(1) ,LIW(1) ,NDG(1) ,NM(1) ,Z(1) 

IF(N.LE.NMAX)GO  TO  10 
WRITE (1,5 )N,NMAX 

5 FORMAT (’  THE  NUMBER  OF  VARIABLES  ',  14 , ’EXCEEDS  THE  STORAGE  LIMIT 
X OF ' , 1 4 , ' N HAS  BEEN  SET  TO  NMAX') 

10  Z(1)=0. 

IF(N.EQ. 1 )GO  TO  25 
DO  20  J=2,N 

IF(NDG(J).GE.NDG(J-1))G0  TO  20 
WRITE (1,15) (NDG (K) ,K=1 ,N ) 

15  FORMAT ( ' THE  VARIABLES  ARE  NOT  ORDERED  ON  DEGREE-PROGRAM  HAS 
X BEEN  TERMINATED* /(IX, 1013)) 

STOP 

20  CONTINUE 
25  1 = 1 

DO  30  J=1 , N 
30  NIW(J)  = I 

DO  55  Msl.MDEG 
NM(M)=I 
LU  =1 
DO  50  J=1,N 
LL  = NIW(J) 

NIW(J)  =1+1 
IF(NDG(J).LT.M)LL  = LL  + 1 
DO  45  LrLL.LU 
IF ( I . LT . KMAX )G0  TO  40 
WRITE (1,35) KM AX 

35  FORMAT ( * THE  NUMBER  OF  COEFFICIENTS  EXCEEDS  THE  STORAGE 
X LIMITS’, 14,'  PROGRAM  HAS  BEEN  TERMINATED') 

STOP 
40  1=1+1 

Z(I )= 10.*Z(L)+J 
45  CONTINUE 
50  CONTINUE 
55  CONTINUE 

NM(MDEG+1 )=I 

RETURN 

END 
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C RMSUB.FOR 

COMPUTATION  OF  U VECTOR;  MODIFIED  DEC  3 1981  FOR  WP 
SUBROUTINE  SETU (N ,MDEG , NIW , LIW, NDG , U , X ) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  NIW( 1 ) , LIW ( 1 ) , NDG (1),U(1),X(1) 
U(1)=1. 

I =1 

DO  5 J=1 ,N 
5 NIW(J)=I 

DO  20  Mrl.MDEG 
LUrl 

DO  15  J = 1,N 
LL  = NIW(J) 

NIW(J)  =1+1 

IF (NDG(J) .LT.M)LL  = LL  + 1 

DO  10  L=LL,LU 

1=1+1 

U(I)=X(J)*U(L) 

10  CONTINUE 
15  CONTINUE 
20  CONTINUE 
RETURN 
END 


CALCULATION  OF  MEANS  AND  SIGMAS,  DATA  NORMALIZATION 

SUBROUTINE  SIGM (N , NDATA , NMAX .MODE , YM , YS , XM , XS , Y , X ) 
IMPLICIT  REAL*8 ( A-H , O-Z ) 

DIMENSION  XM( 1 ) ,XS( 1 ) ,Y( 1 ) ,X(3, 1 ) 

YM  =0 . 

YS=0. 

DO  5 1=1 ,NMAX 
XM(I )=0* 

5 XS(I)=0. 

DO  10  ND=1 , NDATA 
YM=YM+Y(ND) 

YS=YS+Y(ND)*«2 

DO  10  J=1 ,N 

XM  ( J ) =XM  ( J ) +X ( J , ND ) 

10  XS  ( J ) =XS  ( J ) +X  ( J , ND ) * *2 

A=NDATA 
YM=YM/A 
DO  15  J=1 ,N 
XM(J)=XM(J)/A 

15  XS(J)=DSQRT (XS(J)/A-XM( J)**2) 

YS=DSQRT (YS/A-YM##2) 

IF ( MODE . EQ . 0 )G0  TO  25 
DO  20  NDsI, NDATA 
DO  20  J=1 , N 
X ( J , ND ) =X  ( J , ND ) -XM  ( J ) 

20  IF ( MODE ,EQ,2)X(J,ND)=X(J,ND) /XS ( J ) 

25  RETURN 
END 
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C SOLUTION  OF  SYSTEM  EQUATIONS,  SINGLE  ROW  AT  EACH  CALL 
SUBROUTINE  SOLVE(I,A,N,M) 

IMPLICIT  REAL*8 (A-H ,0-Z ) 

DIMENSION  A (20, 1 ) 

IF(I.E0. 1)G0  TO  15 
IM  =1-1 

DO  10  Jrl.IM 

A (I , J ) =0. 

DO  5 K=1 , IM 

5 A(I , J)  = A(I,J)+A(J,K)»A(K,I) 

A(I ,N+1 )=  A(I,N+1)-A(J,I)»A(J,N+1) 

10  A(I,I)  = A(I,I)-A(I,J)*A(J,I) 

15  IF(A(I,I).NE.0.)A(I,I)=1./A(I,I) 

A(I,N+1 )=  A ( I ,N+1 )*A ( I , I ) 

IF(I.EQ.1)G0  TO  30 
DO  25  J = 1 , IM 
A ( J , I ) =-A(I , J )*A(I , I ) 

DO  20  K=J, IM 

A(J,K)  = A(J,K)-A(I,K)*A(J,I) 

20  A(K, J)  = A(J,K) 

A( J ,N+1 )=  A(J,N+1)-A(I,N+1)*A(I,J) 

25  A (I , J ) = A(J,I) 

30  RETURN 

END 
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APPENDIX  C 


LISTING  OF  DELTA  ALGORITHM  FORTRAN  PROGRAM 
FOR  COMPUTATION  OF  DYNAMICALLY  EQUIVALENT  SYSTEM 
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COMPUTATION  OF  DYNAMICALLY  EQUIVALENT  SYSTEM  USING  DELTA  ALGORITHM 
C N JNT=NSEG-1  ONLY 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

BYTE  TITLE(80,6) 

DIMENSION  SEG(22) ,CGS(22) .PHI (3, 3,22) ,CK(3,22) 

DIMENSION  W ( 22 ) , RW  ( 22 ) , DM  ( 22 ) , PH  ( 3 , 22 ) ,BD(6,22) 

DIMENSION  SR (3, 42) ,YPR 1(3,21) , YPR2(3 , 21 ) 

DIMENSION  JOINT (21 ) , JS (21 ) , JNT (21 ) , IPIN (21 ) 

DIMENSION  D(3,3,22),ANG(3,22),PHS(3,3,22) 

C TEMP *********************  DATA  FOR  SAMPLE  CASE 
DATA  LP/2/ 

DATA  DM/-1.25,  0.00,  0.50,-1.00,  0.50,-0.25,  0.25,  0.50, 

X -0.25,  0.25,  0.50,-0.25,  0.50,-0.25,  0.50,7*0.0/ 

C TEMP*********************  END  DATA  FOR  SAMPLE  CASE 
DASIN(Z)  = DTAN2(Z,DSQRT(1 .0D0  - Z*Z)) 

PI  = DTAN2 (0.0D0.-1.0D0) 

RADIAN  = 180.0/PI 

C OPEN  INPUT  FILE,  USES  STANDARD  CVS  FORMAT  FOR  MODEL  DESCRIPTION 
CALL  OPEN (6, 'CVSKANE.DAT  ') 

1 READ (6, 2, END =1 10,ERR=1 10)TITLE 

2 FORMAT (80A1) 

READ(6,5 ) NSEG,NJNT,G 

5 FORMAT(2I6,F10. 0) 

DO  10  I = 1 , NSEG 

10  READ(6 , 15 ) SEG(I ) ,CGS (I ),W(I),(PH(J,I)  ,J=1 ,3) ,(BD(J,I),J=1 ,6) 

15  FORMAT (A4, 1X,A1,10F6.0) 

DO  20  J = 1.NJNT 

20  READ (6, 25 ) JOINT ( J) , JS ( J ) , JNT ( J ) , IPIN ( J ) , (SR ( 1, 2*J-1 ) , 1 = 1 , 3) , 

* (SR(I,2*J),I=1,3).(YPR1(I,J),I=1,3).(YPR2(I,J),I=1,3) 

25  F0RMAT(A4, IX, A1 ,2I4,6F6.0/14X,6F6.0) 

ENTER  WEIGHT  PERTURBATIONS:  (FROM  TERMINAL) 

TEMPORARILY  USE  DATA  STATEMENT  **INSERT  DESIRED  INPUT  STATEMENTS 
WRITE ( 1,26) 

C 26  FORMAT  ( f ENTER  LP,DM(I);  LP  = 2 FOR  PRINTER,  SUM  DM'S  = 07 
C X ’ PROGRAM  WILL  COMPUTE  DM(1)  TO  PRODUCE  A ZERO  SUM’/) 

C READ(1,27)LP,(DM(I), 1=1, NSEG) 

C 27  FORMAT (13, 15F5.0) 

C TEMPORARILY  USE  DATA  STATEMENT  a***************************** 

C END  OF  INPUT  FROM  TERMINAL 
IF(LP. NE. 2)  LP  = 1 
SUM  = O.ODO 
DO  28  I = 2, NSEG 

28  SUM  = SUM  + DM(I) 

DM ( 1 ) = - SUM 
C 

30  F0RMATO0F6.0) 

WRITE  (LP,  31  ) TITLE 

31  FORMAT ( ’ DYNAMICALLY  EQUIVALENT  SYSTEM • /6 ( 1 X , 80A 1 /) ) 

WRITE (LP, 35)  NSEG, N JNT 
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35  FORMATOX,  'SEGMENTS’,  15, ' JOINTS ',15// 

* 32X, ’SEGMENT  MOMENT  OF  INERTIA' , 19X* 'SEGMENT  CONTACT  ELLIPSOID'/ 

* ' SEGMENT ' ,7X, 'WEIGHT ' , 45X, 'SEMIAXES ' , 1 9X, 'CENTER ' /) 

DO  40  I = 1 , NSEG 

40  WRITE (LP, 45 )I,SEG(I),CGS(I) ,W(I ) , (PH ( J , I ) ,J  = 1 ,3) t(BD(J,I) ,J  = 1 ,6) 
45  FORMAT(I3,1X,A4,2X,A1,5X,F7.3,3X,3F11.5,2(5X,3F7.2)) 

WRITE (LP, 50) 

50  FORMAT (/// 

* 3X, 'JOINT' ,15X, 'LOCATION' ,8X, 'SEG(JNT)' , 

* 3X,  'LOCATION ',8X,  'SEGQ+1 )'  , 

* 2X,  'PRIN.  AXIS (DEG)  -SEG(JNT)', 

* 2X,  'PRIN.  AXIS (DEG)  -SEG(J+1)'/ 

* ' J SYM  PLOT  JNT  PIN ' ,2(6X, 'X ' ,8X, 'Y » ,8X, 'Z ' , 3X) , 

» 2(5X, 'YAW' , 5X, 'PITCH' ,5X, 'ROLL* , IX)/) 

DO  55  J = 1 , NJNT 

55  WRITE (LP, 60 )J, JOINT (J) , JS (J) , JNT ( J) , IPIN (J) , (SR (I,2*J-1 ) ,1=1 ,3) , 

* (SR (I , 2*J ) , 1=1 , 3) , ( YPR 1(1, J), 1=1, 3) , ( YPR2(I , J) ,1=1,3) 
60  FORMAT (1X,I2,1X,A4,2X,A1,2I4,4(1X,3F9.2)) 

WRITE(LP,65)(DM(I) ,1=1 ,NSEG) 

65  FORMAT (//'  WEIGHT  PERTURBATIONS' //2(4X, 15F8. 4/)) 

WRITE(LP,66) 

66  FORMAT ( 1H1 ) 

DO  75  N = 1 , NSEG 
DO  75  I = 1,3 
DO  70  J = 1,3 
70  PHKI,  J,N)=0. 

75  PHI (I , I, N)=PH (I , N) 

CALL  DELTA (NSEG, NJNT, JNT, DM, SR, CK, PHI, G,W,RW,LP) 
WRITE(LP,66) 

DO  85  N = 1 , NSEG 
DO  80  I = 1,3 
DO  80  J = 1,3 

80  PHS(I, J,N)=PHI(I, J,N) 

CALL  EIGEN  (PHI ( 1 , 1 ,N ) ,D( 1 , 1 ,N) ) 

ANG(1,N)  = DTAN2(D(2, 1 ,N) ,D(1 , 1 , N))*RADIAN 

ANG(2,N)  = -DAS IN (D( 1 , 3.N ) )*RADIAN 

ANG(3, N)  = DTAN2(D(2,  3,N),D(3,3,N))1‘RADIAN 

WRITE (LP, 90 )N,W(N) , ( ( PHS(I, J,N) , J=1 ,3) ,CK ( I ,N) , PHI (I, I , N) , 

* (D(I,J,N),J=1,3) ,ANG(I , N) , 1=1 , 3) 

85  IF(MOD(N,8).EQ.O)WRm(LP,66) 

90  FORMAT (/'  SEGMENT ', 14 , ' WEIGHT ' ,F7. 3/ 

* 12X, 'INERTIA  MATRIX’ , 12X, 'OFFSET' ,2X, 'EIGENVALUES’ , 6X, 

* 'DIRECTION  COSINE ',9X,'  Y-P-R'/ 

* 3 (5F1 1 .5, 3F10. 6,F 1 1 .5  ^) ) 

WRITE (LP, 95) 

95  FORMAT (1H1, 

* 2X, 'JOINT' ,15X, 'LOCATION' ,8X, 'SEG(JNT)' , 

* 3X, 'LOCATION' ,8X, 'SEG(J+1 )'/ 

* 20X,  2(6X,  'X  ' ,8)(,  'Y  ' ,8X,  'Z  ' , 3X)/) 

DO  100  J = 1.NJNT 
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100  WRITE (LP, 1 05 )J,(SR(I,2*J-1), 1=1,3), (SR(I, 2*J ) , 1=1 , 3) 

105  FORMAT (4X, 1 3, 12X,2( IX, 3F9.2)) 

GO  TO  1 

110  WRITE  (LP, 111) 

111  FORMAT (/'  END  OF  DELTA  PROGRAM.') 

STOP  111 

END 
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SUBROUTINE  DELTA ( NSEG , NJNT , JNT , DM , SR , CK , PHI , G , W , RW , LP ) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  JNT(21),JT(21),IG(22) 

DIMENSION  PHI  ( 3 , 3 . 22 ) , CK ( 3 , 22 ) ,W(22) , RW(22) , DM ( 22* ) , SR  (3,  42 ) 
WRITE(LP, 1) 

1 FORMAT (’  COMPUTATIONS  OF  DELTA  ALGORITHM’/) 

SUM  = 0.0 

DO  5 N = 1 , NSEG 

W(N ) = W(N ) + DM(N) 

RW(N ) = G/W(N ) 

SUM  = SUM  + DM(N) 

DO  5 I =1,3 
5 CK(I , N)=0. 

C DELTA  ALGORITHM 

C 

DO  90  N = 1 , NSEG 

COPY  JOINT  ARRAY 

DO  10  Jsl.NJNT 
10  JT ( J ) = JNT(J) 

SUMT  = 0.0 

DO  35  J = 1 , NJNT 

CHECK  DIRECT  REFERENCE  TO  JOINT 

IF  (JT(J).NE.N)  GO  TO  35 
K = N 
M = 1 

IG(M)  = J + 1 

JM  = 2*J  - 1 

JT (J)  = 0 
SUM  = DM(  J+1 ) 

SUMT  = SUMT  + SUM 

WRITE  (LP,  100)N , J,IC,  DM(  J+1 ) , SUM, SUMT 
100  FORMAT ( IX, 314, 3F10 ,4 ) 

15  DO  20  K=1 , NJNT 

CHECK  FOR  STRING 

IF  (JT(K).NE.IG(M))  GO  TO  20 
M = M + 1 

IG(M)  = K + 1 

SUM  = SUM  + DM(K+1 ) 

SUMT  = SUMT  + DM(K+1) 

WRITE (LP, 100 )N, J, K, DM(K+1 ) , SUM, SUMT 
JT(K)  = 0 i 

20  CONTINUE 

M = M - 1 
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25 
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CHECK  FOR  BRANCHES 


IF  ( M . GT . 0 > GO  TO  15 
»*  PROCESS  JOINT  JM  FOR  RK  ODD  ** 

SRR  = SR  ( 1 , JM) **2  + SR (2 , JM) **2  + SR(3,JM)**2 
DO  30  L = 1,3 

tK(L,N)  = CK(L,N)  + SUM*SR (L , JM) 

DO  25  I = 1,3 

PH I ( I , L , N ) = PHI(I,L,N)  + SUM*SR  (I , JM) *SR (L , JM) /G 
PHI (L, L,N)  = PHI (L, L,N ) - SUM*SRR/G 
CONTINUE 

IF  (N.EQ.1)  GO  TO  50 

ALL  REMAINING  SEGMENTS 

SUM  = -(SUMT+DM(N ) ) 

**  PROCESS  JOINT  N-1  FOR  RK  EVEN  ** 

JM=2*  (N-1 ) 

SRR=SR(  1 , JM)  **2+SR  (2 , JM)  *#2+SR  (3 , JM)  *#2 
DO  45  J = 1,3 

CK( J , N)  =CK(J , N)+SUM*SR ( J , JM) 

DO  40  I = 1,3 

PH I ( I , J , N ) = PHI(I, J,N)  + SUM*SR (I , JM) *SR ( J , JM) /G 
PHI ( J, J,N )=PHI ( J, J,N )-SUM*SRR/G 
CKK  =0. 

DO  55  J = 1,3 
CK(J ,N)=CK(J , N)/W(N ) 

CKK  =CKK+CK  ( J , N ) **2 
CKK=CKK/RW(N) 

DO  60  I = 1,3 
PHI (I , I, N)=PHI (I , I, N)-CKK 
DO  60  J = 1,3 

PH I ( I , J , N ) = PHI(I , J, N)+CK(I , N)*CK(J , N)/RW(N ) 

DO  65  I = 1,3 

IF(N.NE. 1 )SR(I, JM)=SR(I, JM)+CK(I,N) 

DO  75  J = 1 ,N JNT 
IF  ( JNT  ( J ) . NE . N )G0  TO  75 
DO  70  1=1,3 

SR (I , 2*J-1 )=SR(I,2*J-1 )+CK(I , N) 

CONTINUE 

CONTINUE 

RETURN 

END 
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SUBROUTINE  EIGEN  ( A, D) 
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COMPUTES  EIGENVALUES  AND  EIGENVECTORS  OF  A 3X3  SYMMETRIC  MATRIX  A 
EIGENVALUES  ARE  IN  A(I,I),  A(I,J)=0  FOR  I//J. 

EIGENVECTORS  ARE  THE  COLUMNS  OF  D(I,J). 

D IS  THE  DIRECTION  COSINE  MATRIX  RELATING  THE  ORIGINAL  A MATRIX 
TO  THE  PRINCIPAL  SYSTEM  OF  COORDINATES. 

A SUCCESSIVE  ROTATION  ALGORITHM  IS  USED. 

IMPLICIT  REAL*8  ( A-H.O-Z) 

DIMENSION  A(3, 3)  , D(3,3)  , TEMP(3.3) 

TEST  =0.0 
DO  11  1=1,3 

DO  10  J = 1 , 3 

TEST  = TEST  + DABS(A(I,J)) 

D(I, J)  = 0.0 
D(I , I)  = 1.0 
TEST  = TEST*1 . OD-16 
J = 2 
K = 3 
I = 6-K-J 

IF  ( DABS (A(I,J)).LT. TEST ) GO  TO  15 
B = A(I , I)  - A( J , J) 

R = DSQRT(B**2  + 4.0*A(I, J)**2) 

S = DABS(0. 5*B/R) 

C = DSQRTC0.5+S) 

S = DSQRT(0. 5-S ) 

IF  (B*A(I, J) .GT.O.O)  S = -S 
T1  = C*A(I, I)  - S*A(I , J) 

T2  = S*A(I,I)  + C*A(I, J) 

A(I , I)  = C*T 1 - S*(C*A(I, J)-S*A(J, J)) 

A( J, J)  = S*T2  + C*(S*A(I, J)+C*A( J, J) ) 

T1  = C*A(I , K)  - S*A(  J ,K) 

A(J,K)  = S*A(I,K)  + C*A( J,K) 

A(I, K)  = T 1 
DO  14  L = 1 , 3 

T 1 = C#D(L,I)  - 3*D(L,J) 

D(L,  J)  = S *D  ( L , I ) + lC*D(L,  J) 

D(L, I)  = T 1 

A(K, I ) = A(I,K) 

A(K, J)  = A(J,K) 

A(I, J ) = 0.0 

A( J , I)  = 0.0 

IF  ( DABS  ( A ( I , K ) ) +D  A'3S  (A(J,K)).LT.  TEST ) GO  TO  21 
J = 3 
K = K-1 

IF  (K-1)  12,13,13  ! 

DO  23  ITER=1 , 10 
CALL  CFACTT  (D.TEMP, DET) 
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DO  22  1=1,3 

DO  22  J=1 , 3 

D ( I , J ) = 0 . 5 * (D  ( I , J ) +TEMP ( J , I ) /DET ) 

IF  (DABS(D(I, J)) .LT. 1.0D-15)  D(I,J)=0.0 

IF  (DABS (DET- 1 .0) . LT. 1 . OD-6 ) GO  TO  25 

23  CONTINUE 

WRITE  (1,24)  DET  i 

24  FORM AT ( '0  EIGEN  RENORMALIZATION  DID  NOT  CONVERGE,  DET 
X , 1PD25. 15 ) 

25  RETURN 
END 
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SUBROUTINE  CFACTT  ( A,B,D) 


GIVEN  3X3  MATRIX  A 

COMPUTE  B TRANSPOSE  OF  COFACTORS  (SIGNED  MINORS) 
AND  D THE  VALUE  OF  THE  DETERMINANT  OF  A. 

INVERSE  OF  A IS  B(J,K)/D 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  A( 3, 3) ,B(3,3) 

M = 4 
L = 2 
N r 3 
D = 0.0 
DO  20  J=1 , 3 

B( J, J)  = A(L,L)*A(NfN)-A(L,N)*A(N,L) 

IF  (J.EQ.3)  GO  TO  20 
L = N 
N = J 
KK  = J+1 
DO  15  K=KK, 3 
M = M-1 

B(K, J)  = A(K,M)*A(M,J)-A(K, J)*A(M,M) 

B(J,K)  = A(J,M)*A(M,K)-A(J,K)*A(M,M) 

D = D+A( 1 , J)*B( J , 1 ) 

RETURN 

END 


